O kolik déle čekáme po delší erupci — a jak jistě to víme?

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.

Logo kurzu Biostatistika BIOSTATISTIKA

L04

Od odhadu k nejistotě

cuni-natur-biostatistics.github.io/materialy.html

Ilustrace vytvořená pomocí AI.

Gejzír Old Faithful v číslech

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.

Léto 2024

  • 117 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


Napoví délka právě skončené erupce, jak dlouho budeme čekat na další?

Foto: NPS/Jim Peaco, public domain. Rozsahy: National Park Service.

Co si pamatujete z minulé lekce?

Odpovězte samostatně na tři krátké otázky. Společné odpovědi projdeme po hlasování.

Dva kosatce se liší šířkou korunního lístku o 0,5 cm. O kolik se podle modelu v průměru liší jejich odhadnutá délka lístku?

Výstup R příkazu coef(mod_listky): intercept 1.08 a sklon sirka_listku 2.23.

  • Asi o 2,2 cm
  • Asi o 4,5 cm
  • Asi o 1,1 cm
  • Nelze určit bez šířek obou květů

Správná odpověď: Asi o 1,1 cm

Sklon 2,23 cm/cm vynásobíme rozdílem šířek: 2,23 · 0,5 ≈ 1,1 cm. Intercept se při porovnání dvou květů odečte, takže jejich konkrétní šířky nepotřebujeme.

Residua kosatců druhu versicolor leží převážně nad nulou. Co z toho plyne?

Graf residuí proti odhadnuté délce lístku, body barevně podle druhu: setosa leží převážně pod nulou, versicolor převážně nad nulou a virginica na obou stranách.

  • Model jejich délku převážně podhodnocuje; v modelu chybí informace o druhu.
  • Model jejich délku převážně nadhodnocuje; v modelu chybí informace o druhu.
  • Nic; residua všech květů mají dohromady průměr nula, model je v pořádku.
  • Versicolor má ze všech druhů nejdelší lístky, proto má kladná residua.

Správná odpověď: Model jejich délku převážně podhodnocuje; v modelu chybí informace o druhu.

Kladné residuum: naměřená délka je větší než odhad, model ji podhodnocuje. Průměr nula platí pro všechna residua dohromady, ne pro každý druh; posun skupiny ukazuje chybějící prediktor.

Dvě přímky mají pro stejné čtyři květy tato residua. Která přímka popisuje květy lépe?

Tabulka residuí čtyř květů v centimetrech: první přímka +0,5, −0,5, +0,5, −0,5; druhá přímka +0,2, −0,1, +0,2, −0,1.

  • První přímka, protože součet jejích residuí je nula.
  • Obě stejně, protože obě mají kladná i záporná residua.
  • Nelze rozhodnout, dokud neznáme intercept a sklon obou přímek.
  • Druhá přímka, protože součet čtverců jejích residuí je menší.

Správná odpověď: Druhá přímka, protože součet čtverců jejích residuí je menší.

Kladná a záporná residua se v součtu ruší, proto je umocňujeme: první přímka má součet čtverců 1,00 cm², druhá jen 0,10 cm².

Výsledky učení


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

  • vysvětlit, proč jiný soubor pozorování dává jiný odhad sklonu
  • najít ve výstupu summary() odhad sklonu a jeho standardní chybu, vysvětlit, co standardní chybu zmenšuje, a rozlišit velikost vztahu od přesnosti odhadu
  • získat 95% interval spolehlivosti sklonu funkcí confint() a pás spolehlivosti 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, zda jedna erupce výrazně mění odhad sklonu
  • zapsat výsledek jednou větou s odhadem, SE, intervalem a jednotkami

Dají jiné erupce stejnou přímku?

Každý řádek je jedna erupce

# 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)
[1] 117
# Celkový počet chybějících hodnot
sum(is.na(x = data_gejzir))
[1] 0

old_faithful_2024.csv: uložte do složky data ve svém projektu v R.

Delší erupce, delší čekání?

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í.

Jaký směr má vztah? Kde je nejvíc bodů? Leží body na jediné přímce?

Přímka z tohoto léta

Stejný bodový graf s fialovou rostoucí přímkou lineárního modelu.

Odhadnutý sklon: 13,34 minuty čekání na minutu erupce

Stejná data + stejný vzorec modelu = vždy stejné odhady. Funkce lm() nic náhodně nelosuje.

Sedí přímka k datům?

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.

Histogram residuí: většina residuí leží blízko nuly, několik je výrazně kladných.

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.

Každý model zkontrolujte, než ho začnete vykládat, i ve vlastních datech: standardní chyba a interval platí jen pro přímku, která data přiměřeně popisuje.

Dostali bychom příště stejnou přímku?

Otázka

Gejzír pozoruje někdo jiný v jiných dnech, zaznamená své body a fituje stejný model. Jaký sklon dostane?


A) Přesně stejný sklon jako my.

B) Podobný, ale ne stejný sklon.

C) Úplně jiný sklon, nesouvisející s naším.


Zvolte A, B, nebo C a zdůvodněte svou volbu sousedovi.

Pepa pozoroval pouze červen

Bodový graf Pepových červnových erupcí s fialovou přímkou odhadnutou pouze z nich.

V Pepových datech se dvě erupce, jejichž délky se liší o minutu, v průměru liší čekáním o 12,92 minuty.

Co čekáte od Mařenky a Karla?

Diskuze se sousedem — 1 min

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?

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.

Ilustrace vytvořená pomocí AI.

Stejný postup, tři různé výsledky

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í.

Rozdíly mezi měsíci mohou odrážet jiné zaznamenané erupce i skutečnou změnu chování gejzíru. Tři měsíce tyto příčiny neoddělí.

Odhad závisí na datech


Co zůstalo stejné

  • otázka a proměnné
  • lineární model
  • funkce lm()

Co se změnilo

  • pozorované erupce
  • odhadnutá přímka a její sklon


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 (sampling variability).

Sto pozorovatelů, sto sklonů

Skutečnost, výběr, odhad

Ve výzkumu

Skutečnost

Jak se gejzír chová u všech erupcí (cílová populace). Skutečný sklon nevidíme.

→

Výběr

Erupce, které někdo zaznamenal: našich 117.

→

Odhad

Jeden sklon: 13,34.

V simulaci

Skutečnost si zvolíme

Přímka z tohoto léta jako referenční sklon.

→

100 výběrů

Každý ze 100 pozorovatelů: 40 erupcí.

→

100 odhadů

Sto sklonů.

V simulaci víme, kde je pravda; ve výzkumu máme jeden výběr a skutečnost neznáme.

Co udělá sto dalších pozorovatelů?

Každý pošle do aplikace GeyserWatch 40 erupcí, tedy jeden výběr, odhadne stejný model a získá jeden sklon.

Jak budou vypadat všechny sklony společně?

Ve dvojici nakreslete hrubý tvar, který očekáváte po 100 hlášeních.

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.

Ilustrace vytvořená pomocí AI.

Simulace

Délky erupcí z našich dat; čekání vytvoříme z referenční přímky (sklon 13,34) s odchylkami velkými jako residua.

Sledujte cestu jednoho sklonu

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.

Co se mění z hlášení na hlášení a co zůstává stejné?

Nový výběr bodů → nová přímka → nový sklon → další díl histogramu.

Po 100 hlášeních

Poslední snímek animace: body posledního hlášení se všemi předchozími přímkami v pozadí a úplný histogram sta sklonů.

Co vidíme — a co ještě nevíme?


Co vidíme

  • sklony nejsou stejné
  • většina leží ve společné oblasti
  • některé výběry dávají strmější, jiné plošší přímku

Co ještě nevíme

  • jak šířku histogramu 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í?

Histogram sta simulovaných sklonů se svislou tmavou čarou s popiskem referenční sklon.

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.

Ze sta hlášení

Směrodatná odchylka 100 sklonů: 1,095

Z jediného hlášení

R vypíše standardní chybu (SE), typicky 1,091

Standardní chyba (SE) odhaduje z jediného výběru, jak moc by se sklony typicky lišily mezi výběry stejné velikosti.

Co dělá sklon stabilnějším?

Otázka

100 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?


Srovnání 1 · body kolem přímky

A) body daleko od přímky

B) body blízko přímky

Srovnání 2 · délky erupcí

A) délky z celého rozsahu

B) jen erupce od 3 minut

Srovnání 3 · počet erupcí

A) 10 erupcí v hlášení

B) 40 erupcí v hlášení


Se sousedem zvolte v každém srovnání A, nebo B. Na odpověď máte minutu.

Tři srovnání, tři odpovědi

Š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ší.

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é.

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é.

Stabilnější sklon: body blízko přímky, délky z celého rozsahu (body daleko od sebe drží přímku jako houpačku), více erupcí. Malý histogram v rohu: užší znamená menší SD sklonů.

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}} \]


1 · Čitatel

Body daleko od přímky → větší SE

2 · Jmenovatel

Délky z celého rozsahu → menší SE

3 · Jmenovatel

Více erupcí → menší SE, ale jen s odmocninou

Jak velké je typické residuum?

Residua jsou svislé odchylky erupcí od přímky. Jejich typickou velikost shrneme podobně jako SD.


1 · Sečteme čtverce

RSS = 5368,7 min²

Součet čtverců residuí (v minulé lekci SSE): lm() volí přímku s nejmenším RSS.

2 · Vydělíme n − 2

5368,7 / 115 = 46,68 min²

3 · Odmocníme

√46,68 = 6,83 min

\[ \text{residuální směrodatná odchylka} = \sqrt{\frac{\text{RSS}}{n - 2}} = \sqrt{\frac{\text{5368,7}}{115}} = \text{6,83}\ \text{min} \]

Residuální směrodatná odchylka (ve výstupu R Residual standard error): erupce leží typicky asi 6,83 min nad přímkou nebo pod ní.

Spočítejte standardní chybu sami

Všech 117 erupcí

  • residuální směrodatná odchylka: 6,83 min
  • SD délek erupcí: 0,965 min
  • počet erupcí: 117

Dvojice — 2 minuty

\[ \text{SE} = \frac{\text{residuální SD}}{\text{SD délek} \times \sqrt{n - 1}} \]

Kalkulačka v telefonu stačí.

\[ \begin{aligned} \text{SE} &= \frac{\text{6,83}} {\text{0,965} \times \sqrt{117 - 1}} \\[4pt] &= \frac{\text{6,83}}{\text{10,39}} = \text{0,66} \end{aligned} \]

SE sklonu = 0,66 min/min — z jediného souboru dat.

Proč má hlášení se 40 erupcemi větší SE?

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.

Stejná křivka s dalším bodem pro čtyřnásobek našich erupcí, kde je standardní chyba poloviční.

Hlášení: vzorec 1,13, simulace 1,09. Čtyřikrát víc erupcí → jen poloviční SE.

Co vypíše summary()?

# Model ze všech erupcí a jeho souhrn
mod_gejzir <-
  lm(
    formula = cekani_min ~ delka_erupce_min,
    data = data_gejzir
  )

summary(object = mod_gejzir)

Call:
lm(formula = cekani_min ~ delka_erupce_min, data = data_gejzir)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.0280  -4.4756  -0.3648   3.6354  21.0827 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       45.9122     2.4658   18.62   <2e-16 ***
delka_erupce_min  13.3350     0.6575   20.28   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.833 on 115 degrees of freedom
Multiple R-squared:  0.7815,    Adjusted R-squared:  0.7796 
F-statistic: 411.4 on 1 and 115 DF,  p-value: < 2.2e-16

Dvojice — 1 min

Co z výstupu už znáte? Co vidíte poprvé?

  • 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
  • zbytek: v dalších lekcích

Tabulka Coefficients: velikost a přesnost

# Jen tabulka Coefficients ze souhrnu modelu, na 3 platné číslice
signif(
  x = summary(object = mod_gejzir)$coefficients,
  digits = 3
)
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)          45.9      2.470    18.6 1.66e-36
delka_erupce_min     13.3      0.657    20.3 8.67e-40


Dvojice — 1 min

Najděte číslo, které popisuje

  1. velikost vztahu,
  2. přesnost tohoto odhadu.

Řádek delka_erupce_min

  • Estimate = 13,34 min/min: velikost vztahu
  • Std. Error = 0,66 min/min: stejné číslo, jaké jste spočítali

Velikost vztahu, nebo přesnost?

Ve dvojici přiřaďte každou otázku k jednomu ze dvou čísel.

„O kolik minut se v průměru liší čekání po erupcích lišících se o minutu?“

Velikost vztahu (Estimate)

„Jak moc by se odhad lišil mezi výběry?“

Přesnost odhadu (Std. Error)

„Jak těsně se možné odhady shlukují?“

Přesnost odhadu (Std. Error)

„Je vztah prudký, nebo pozvolný?“

Velikost vztahu (Estimate)

O kolik déle čekáme — a jak jistě to víme?


Už víme

  • velikost vztahu: o 13,34 min čekání na minutu erupce
  • přesnost odhadu: SE = 0,66 min/min
  • SE klesá s více erupcemi a s širším rozsahem délek

Ještě nevíme

  • jaký rozsah sklonů je s daty slučitelný
  • jak přesně známe celou přímku
  • zda závěr nestojí na jediné erupci

Jaké sklony jsou s daty slučitelné?

Otázka

Každé ze 100 hlášení GeyserWatch si kolem svého sklonu vyznačí rozsah odhad ± asi 2 SE. Při mnoha opakováních zachytí referenční sklon kolik takových intervalů?


A) asi 50

B) asi 70

C) asi 95

D) všech 100


Zvolte A, B, C, nebo D a zdůvodněte svou volbu sousedovi.

Sto pozorovatelů, sto 95% intervalů

Stejných 100 hlášení, každé se svým rozsahem odhad ± asi 2 SE (přesně: 95% interval z confint()).

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í.

97 ze 100 95% intervalů zachytilo referenční sklon. Ve skutečném výzkumu skutečný sklon neznáme.

95 % patří k postupu

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.

Ilustrace vytvořená pomocí AI.

95% interval spolehlivosti (confidence interval, CI)

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).

Neříkáme

  • 95 % pozorování je uvnitř
  • tento interval má 95% šanci obsahovat pravdu
  • interval dokazuje určitou hodnotu

Jak takový interval vytvořit z jednoho výběru?

Máme odhad 13,335 a SE 0,657 (na tři desetinná místa). Kolik standardních chyb přidat na každou stranu?

\[ \text{polovina šířky} = \text{kritická hodnota} \times \text{SE} = \text{1,981} \times \text{0,657} = \text{1,302} \]

\[ \text{dolní hranice} = \text{13,335} - \text{1,302} = \text{12,033} \]

\[ \text{horní hranice} = \text{13,335} + \text{1,302} = \text{14,637} \]

Pro 95% interval je kritická hodnota blízko 2; její přesnou hodnotu vysvětlíme v následující lekci.

Stejný interval v R


# 95% intervaly spolehlivosti koeficientů
confint(object = mod_gejzir)
                    2.5 %   97.5 %
(Intercept)      41.02795 50.79642
delka_erupce_min 12.03274 14.63731

Řádek delka_erupce_min

12,03 až 14,64 minuty čekání na minutu erupce — stejné hranice jako ručně.


Sloupce 2.5 % a 97.5 %: 95% interval nechává na každé straně 2,5 %.

95% interval spolehlivosti

Vodorovný fialový interval sklonu s tmavým bodem odhadu uprostřed a popisky dolní hranice, odhadu a horní hranice.

Pro Old Faithful: [12,03; 14,64] minuty čekání na minutu erupce.

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


Dvojice — 90 sekund

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.


Jedna možná oprava

„Při mnoha opakováních zachytí takto vytvořené intervaly skutečný sklon přibližně v 95 % případů.“

Co se stane s 80% intervaly?

Otázka

Stejných 100 hlášení, ale místo 95% intervalů spočítáme 80%. Jaké budou intervaly a kolik jich zachytí skutečný sklon?


A) Užší, a zachytí ho méně z nich.

B) Širší, a zachytí ho více z nich.


Ve dvojici zvolte A, nebo B a připravte jedno vysvětlení.

Sto pozorovatelů, nižší hladina spolehlivosti

Stejných sto hlášení s osmdesátiprocentními intervaly: intervaly jsou kratší a více z nich je šedých.

Stejných sto hlášení s padesátiprocentními intervaly: intervaly jsou krátké a zhruba polovina z nich referenční čáru míjí.

80 %: 78 ze 100 intervalů zachytilo referenční sklon. 50 %: 54 ze 100. Hladina spolehlivosti popisuje úspěšnost postupu.

Užší neznamená automaticky přesnější

Čtyři vodorovné intervaly pod sebou se stejným středem; osmdesátiprocentní je nejkratší, devadesátidevítiprocentní nejdelší.

Intervaly porovnáváme jako přesnější a méně přesné jen při stejné hladině spolehlivosti.

Snížením hladiny z 95 % na 80 % získáme užší interval, ale z dat jsme nezískali víc informace: data i SE zůstaly stejné.

Jak přesně známe celou přímku?

Otázka

Průměrné čekání pro každou délku erupce, celý tento vztah, nazýváme průměrný trend; přímka ho odhaduje. Nakreslíme přímky všech 100 hlášení GeyserWatch do jednoho grafu. Kde se od sebe budou lišit nejméně?


A) U krátkých erupcí.

B) U průměrně dlouhých.

C) U dlouhých erupcí.


Zvolte A, B, nebo C a zdůvodněte svou volbu sousedovi.

Sto hlášení, sto přímek

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íř.

Fialové jsou přímky hlášení, tmavá je přímka ze všech erupcí. Přímky se kříží u průměrné délky erupce (čárkovaně) a k okrajům se rozbíhají.

Pás pro průměrný trend

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žší.

Pás spolehlivosti: 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ů.

Kdy je pás široký?

Č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.

Pás rozšiřují stejné ingredience jako SE: velký rozptyl kolem přímky, málo bodů, úzký rozsah x. Nejširší je daleko od dat: tam přímku prodlužujeme naslepo.

Kolik erupcí leží uvnitř pásu?

Bodový graf všech erupcí s fialovou přímkou a světle fialovým pásem spolehlivosti.

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.

Pás ≠ kde leží 95 % erupcí

Uvnitř leží jen 22 ze 117 erupcí. Pás ukazuje, kde může ležet průměrný trend. Rozsah pro jedno nové čekání musí být širší.

Může jedna erupce změnit závěr?

Otázka

Při přípravě dat jsme tři erupce vyřadili. Která by po vrácení naklonila přímku nejvíc?


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.

A) 2,02 min erupce, čekání 94 min

B) 4,27 min erupce, čekání 121 min

C) 0,68 min erupce, čekání 106 min

Krajní erupce mají delší páku

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.

Pákový efekt (leverage)

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í.

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

Krajní erupce je jako dítě na konci houpačky: přímku nakloní nejsnáz.

Kdy jeden bod pohne přímkou?

Č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í.

Cookova vzdálenost: jak moc by se změnily odhadnuté hodnoty bez bodu. Velká je až při velké páce a velkém residuu současně.

Delší páka a velké residuum, velký vliv

Cookova vzdálenost (Cook’s distance)

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í.

# Cookova vzdálenost každé erupce
cooks.distance(model = mod_gejzir_120)
Erupce Residuum (min) Pákový efekt Cookova vzdálenost
A 18,2 0,029 0,072
B 18,4 0,012 0,029
C 46,1 0,079 1,387

A a B jsou odlehlé se stejně velkým residuem; 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.

Mění některá z našich erupcí závěr?

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á.

Čára 4/n

4 děleno počtem erupcí: orientační hodnota, ne hranice pro mazání. V modelu ze 117 erupcí ji přesahuje 3 z 117 erupcí: na ty se podíváme blíž.

Vynecháme jednu erupci

Sklon se změní nejvýš o 0,18 min/min, tedy asi o 27 % SE (0,66).

Závěr nestojí na jediné erupci: vliv jedné erupce porovnáváme se standardní chybou.

Co s vlivným pozorováním?

Prověřit, ne mazat

  • 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í.

Naše tři erupce

  • A a B posunou sklon jen o 0,20 min/min (z 13,34 na 13,14). 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 11,92 (asi 1,9 × SE našeho modelu). Trvala jen 41 s, tedy pod rozsahem 1,5–5 min podle National Park Service (NPS): důvod k vyřazení máme.

Závěr

Co už víme o čekání na další erupci?

O kolik déle čekáme po delší erupci — a jak jistě to víme?


Co už můžeme říct

Dvě erupce, jejichž délky se liší o minutu, se v průměru liší čekáním o 13,34 min; s daty a modelem jsou slučitelné sklony od 12,03 do 14,64.

Po čtyřminutové erupci je odhadnuté průměrné čekání 99,25 min (95% CI [97,91; 100,60]).

Co zůstává otevřené

Jak pomocí dat posoudit konkrétní tvrzení o sklonu — například že délka erupce s čekáním vůbec nesouvisí?

Napište jednu vědeckou větu

Dvojice — 3 min

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ě

Data a výsledek modelu

  • 117 erupcí gejzíru Old Faithful, léto 2024
  • délka erupce: 1,65–4,83 min
  • čekání: 62–116 min
  • odhad sklonu: 13,34
  • SE: 0,66
  • 95% CI: [12,03; 14,64]
  • jednotky: min čekání / min erupce

Vyměňte si větu se sousední dvojicí: jsou v ní proměnné, jednotky a nejistota?

Jedna možná věta

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 13,34 minuty (odhad sklonu = 13,34, SE = 0,66, 95% CI [12,03; 14,64] minuty čekání na minutu erupce).

Co ještě omezuje náš závěr?

Dobrovolníci nezaznamenali náhodný výběr všech erupcí a blízké erupce mohou být časově propojené.


Výsledek je opatrný odhad vztahu v zachyceném období — ne neměnný přírodní zákon. Nestojí ale na jediné erupci.

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.

Ilustrace vytvořená pomocí AI.

Shrnutí

  • Jiný výběr dává jiný odhad sklonu: výběrová variabilita.
  • Std. Error je přesnost odhadu: roste s residui, klesá s rozsahem prediktoru a s odmocninou z počtu pozorování.
  • 95% interval a pás spolehlivosti: 95 % patří k postupu, ne k pozorováním.
  • Vlivné pozorování = dlouhá páka a velké residuum; prověřit, ne mazat.
  • Vědecká věta: odhad, SE, interval a jednotky.

Co potřebujete vidět, než napíšete „po minutě erupce navíc čekáme o 13,34 minuty déle“?


Jaké tři otázky položíte autorovi, který vám ukáže jen tento titulek?


Řekněte sousedovi jednou větou

Kterou z otázek byste položili jako první, a proč?

Jak posoudit konkrétní tvrzení o sklonu? Následující lekce přidá testování hypotéz.