Jakou délku korunního lístku odhadneme z jeho šířky?

Lineární model s jedním numerickým prediktorem

1 Úvod

V minulé lekci jsme popisovali vztah délky ploutve a hmotnosti tučňáků pomocí bodového grafu a korelace. Korelace řekla, jak těsně dvě měření souvisejí, ale neodpověděla na otázku „o kolik“. V této lekci tuto otázku položíme korunním lístkům kosatců: obě měření ponecháme v původních centimetrech a pomocí lineárního modelu odhadneme délku lístku z jeho šířky.

QuestionHlavní otázka této lekce

Mají širší korunní lístky kosatců také větší délku, a pokud ano, jak velká je tato změna?

Biologická otázka určuje role měření. Korelace se po prohození os nezmění; model pro odhad délky z šířky však řeší jinou otázku než model pro odhad šířky z délky. Délka korunního lístku je proto odezva a šířka prediktor.

QuestionCo čekáte ještě před prvním grafem?

• A: Širší lístky budou téměř vždy delší, a to u všech květů o stejný kus.
• B: Širší lístky budou v průměru delší, ale jednotlivé květy se budou od tohoto průměrného vztahu lišit.
• C: Šířka a délka korunního lístku spolu souviset nebudou.

Než budete číst dál, vyberte možnost, která vám připadá nejpravděpodobnější, a zdůvodněte ji.

DůležitéOdpověď: očekáváme průměrný, ale nedokonalý vztah

Biologicky nejpravděpodobnější je možnost B. Rozměry částí jednoho květu spolu obvykle souvisejí, ale květy se liší druhem, věkem i podmínkami, ve kterých rostly. Model proto popíše průměrný vztah a u jednotlivých květů zůstanou odchylky.

2 Výsledky učení

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

• převést biologickou otázku na vztah mezi odezvou a prediktorem
• fitovat model lm(odezva ~ prediktor) a interpretovat směr i velikost efektu
• vysvětlit residuum jako rozdíl mezi naměřenou a odhadnutou hodnotou
• číst graf residuí a rozpoznat v něm zbývající vzor

3 Každý květ má dvě měření

Pro tuto lekci používáme reálná měření kosatců, která sesbíral botanik Edgar Anderson a později je publikoval Ronald Fisher. Jeden řádek tabulky odpovídá jednomu květu. Šířka a délka korunního lístku jsou numerické proměnné v centimetrech; druh kosatce je kategoriální proměnná.

V prvním modelu použijeme záměrně jen jeden numerický prediktor a druh zatím do modelu nepřidáme. Na konci se k tomuto rozhodnutí vrátíme při diagnostice.

TipDoplňující: proč kosatce?

• měří biologicky srozumitelné veličiny ve stejných jednotkách
• obsahují přirozenou variabilitu mezi květy
• hodnoty šířky začínají blízko nuly, takže intercept můžeme ukázat bez stlačení dat do malé části grafu
• umožňují plynulý přechod od bodového grafu k lineárnímu modelu

Soubor kosatce.csv je součástí materiálů. Pro následující kód jej uložte do složky data ve svém projektu v R. Načtením vznikne tabulka data_kosatce; funkce str() ukáže názvy, typy a první hodnoty sloupců:

# Načteme tabulku kosatců ze souboru CSV
data_kosatce <-
  read.csv(file = "data/kosatce.csv")

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

'data.frame':   150 obs. of  4 variables:
 $ sirka_listku          : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ delka_listku          : num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ druh                  : chr  "setosa" "setosa" "setosa" "setosa" ...
 $ delka_kalisniho_listku: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...

Tabulka má čtyři sloupce:

• sirka_listku: šířka korunního lístku v centimetrech (prediktor)
• delka_listku: délka korunního lístku v centimetrech (odezva)
• druh: druh kosatce, který zatím ponecháme mimo model
• delka_kalisniho_listku: délka kališního lístku v centimetrech; v této lekci ji nepoužijeme

Dále zkontrolujeme rozsahy hodnot, chybějící hodnoty a počet květů jednotlivých druhů:

# Minimum, kvartily, průměr a maximum numerických sloupců
summary(
  object = data_kosatce[, c("sirka_listku", "delka_listku")]
)

  sirka_listku    delka_listku  
 Min.   :0.100   Min.   :1.000  
 1st Qu.:0.300   1st Qu.:1.600  
 Median :1.300   Median :4.350  
 Mean   :1.199   Mean   :3.758  
 3rd Qu.:1.800   3rd Qu.:5.100  
 Max.   :2.500   Max.   :6.900  

# Počet chybějících hodnot v každém sloupci
colSums(x = is.na(x = data_kosatce))

          sirka_listku           delka_listku                   druh 
                     0                      0                      0 
delka_kalisniho_listku 
                     0 

# Počet květů jednotlivých druhů
table(data_kosatce$druh)

    setosa versicolor  virginica 
        50         50         50 

V datech je 150 květů ze 3 druhů. Šířky se pohybují od 0,1 do 2,5 cm a délky od 1 do 6,9 cm; žádný sloupec nemá chybějící hodnoty. Jsou to malé kladné rozměry odpovídající květům, takže nic nenaznačuje chybu jednotek ani nesmyslnou hodnotu.

TipDoplňující: jak vznikl soubor kosatce.csv

Měření korunních a kališních lístků tří druhů kosatců (Iris setosa, I. versicolor a I. virginica) sesbíral Edgar Anderson (Anderson 1935, The irises of the Gaspé Peninsula, Bulletin of the American Iris Society 59: 2–5). Ronald Fisher je použil v článku o rozlišování druhů pomocí několika měření (Fisher 1936). Data jsou součástí základní instalace R jako tabulka iris v balíčku datasets; podrobnosti uvádí dokumentace datasetu.

Soubor pro tuto lekci obsahuje všech 150 řádků tabulky iris. Žádný řádek ani naměřenou hodnotu jsme nezměnili; vybrali jsme jen čtyři sloupce a dali jim české názvy: sirka_listku (původně Petal.Width), delka_listku (Petal.Length), druh (Species) a delka_kalisniho_listku (Sepal.Length). Protože je R šířeno pod licencí GNU GPL verze 2 nebo 3, šíříme za stejných podmínek i tento výběr.

Soubor vytvořil skript R/prepare_l03_data.R, který před uložením ověří počet řádků a absenci chybějících hodnot. Stejnou tabulku si můžete vytvořit i sami v základním R:

data_export <-
  data.frame(
    sirka_listku = iris$Petal.Width,
    delka_listku = iris$Petal.Length,
    druh = as.character(iris$Species),
    delka_kalisniho_listku = iris$Sepal.Length
  )
write.csv(
  x = data_export,
  file = "data/kosatce.csv",
  row.names = FALSE
)

Pro výuku používáme samostatný soubor CSV, protože v praktiku stačí načíst jediný soubor pomocí read.csv().

3.1 Šířka a délka jednotlivých květů

Než budeme mluvit o modelu, podívejme se na jednotlivé květy.

Z předchozích lekcí už víme dvě základní věci:

• obě měřené proměnné jsou numerické, tedy měřené na číselné škále
• každou proměnnou můžeme shrnout samostatně, ještě než se podíváme na jejich vztah

Porovnejte rozsahy obou měření. Stačí jejich samostatná shrnutí k určení toho, zda má širší lístek také větší délku? Doplňující grafy ukazují každou proměnnou samostatně; vztah vyžaduje zachování dvojic ze stejného květu.

TipDoplňující: zopakování každé proměnné samostatně

Šířka korunních lístků

Jaké šířky jsou pro tyto květy běžné a jaké krajní? Porovnejte minimum, medián a maximum se zobrazenými měřeními.

# Zkontrolujeme polohu a rozptýlení šířek korunních lístků
summary(object = data_kosatce$sirka_listku)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.100   0.300   1.300   1.199   1.800   2.500 

Délka korunních lístků

Mají délky stejný rozsah jako šířky? Obě veličiny měříme v cm; z rozdílných rozsahů samotných ještě nepoznáme vztah v jednotlivých květech.

# Zkontrolujeme polohu a rozptýlení délek korunních lístků
summary(object = data_kosatce$delka_listku)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.600   4.350   3.758   5.100   6.900 

Jaký typ grafu dává největší smysl, když máme dvě numerické proměnné měřené na stejném květu?

DůležitéOdpověď: pro dvě měření zachováme jejich dvojice

Dobrá volba je bodový graf, protože chceme vidět jednotlivé květy a vztah mezi dvěma měřeními.

3.2 Jak spolu souvisejí šířka a délka?

Prohlédněte si body: roste délka s šířkou, je průběh přibližně přímý a leží všechny květy stejně blízko trendu? Svůj popis porovnejte s následujícími vodítky.

V našem případě je prediktorem šířka korunního lístku a podle ní vysvětlujeme rozdíly v délce korunního lístku.

Co z grafu sledujte zejména:

• směr vztahu: zda se s rostoucí šířkou obvykle zvětšuje i délka
• neobvyklé květy: jestli jsou v datech květy, které se chovají výrazně jinak než ostatní
• rozptýlení kolem trendu: jak daleko leží jednotlivé hodnoty od hlavního trendu

Jednotlivé květy v

plot(
  # Každý květ vytvoří jeden bod se dvěma naměřenými souřadnicemi
  x = data_kosatce$sirka_listku,
  y = data_kosatce$delka_listku,
  # Popisky vysvětlují význam a jednotky obou os
  xlab = "Šířka korunního lístku (cm)",
  ylab = "Délka korunního lístku (cm)",
  main = "Bodový graf šířky a délky korunního lístku",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.4,
  cex.axis = 1.5,
  cex.lab = 1.5,
  cex.main = 1.5
)

Graf pomocí {ggplot2}

ggplot2::ggplot(
  # Data určují řádky; aes() přiřadí sloupce k osám
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  # Jeden řádek dat se zobrazí jako jeden bod
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  # Popíšeme osy i celý graf
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)",
    title = "Bodový graf šířky a délky korunního lístku"
  ) +
  ggplot2::theme_bw(base_size = 18)

3.3 Která přímka vystihuje měření?

Vyberte z grafu jednu přímku a zdůvodněte svou volbu. Které vzdálenosti bodů od přímky by mělo pravidlo pro výběr zohlednit?

Pro stejná data existuje více kandidátních přímek a každá vede k trochu jiné interpretaci. Každá přímka má přitom jiné parametry (viz níže).

Než rozhodneme, která přímka vystihuje data nejlépe, potřebujeme umět libovolnou přímku jednoznačně popsat.

4 Jak popíšeme polohu přímky?

Každou přímku popíšeme dvěma parametry:

• a: intercept (anglicky intercept; česky také průsečík)
• b: sklon (anglicky slope; česky také regresní koeficient)

Intuitivně:

• a posouvá přímku nahoru nebo dolů (mění základní úroveň odhadnuté délky)
• b mění strmost přímky, tedy jak rychle roste odhadnutá délka s rostoucí šířkou

Před prohlédnutím grafů odhadněte: který parametr mění výšku přímky a který rozdíl mezi jejími odhady při dvou šířkách? V prvním grafu se mění pouze intercept, v druhém pouze sklon.

Co se změní, když přímku posuneme?

Co se změní, když přímku nakloníme?

Obrázek níže ukazuje čtyři kombinace obou parametrů nad stejnými květy. V řádcích se mění intercept a, ve sloupcích sklon b.

Každou kandidátní přímku tedy popíšou dvě čísla. Jak ale poznáme, která z nich nejlépe vystihuje vztah mezi šířkou a délkou korunního lístku?

PoznámkaJak daleko leží naměřené délky od přímky?

Potřebujeme kritérium, které využije vzdálenosti naměřených délek od přímky při dané šířce. Proto potřebujeme odhadnuté délky a residua.

5 Čtyři květy a jedna zvolená přímka

5.1 Čtyři květy napříč prostorem dat

Než porovnáme kandidátní přímky ve všech datech, ukážeme princip na čtyřech květech rozprostřených napříč grafem.

TipČtyři skutečné květy pro ruční výpočet

Vybrané řádky zachovávají dvojice měření stejného květu. Čísla řádků určují příklad; nejsou součástí statistické metody.

# Vybereme čtyři květy a ponecháme jen dvě modelované proměnné
data_ctyri_kvety <-
  data_kosatce[
    c(3, 82, 117, 101),
    c("sirka_listku", "delka_listku")
  ]

# Zobrazíme malou tabulku, se kterou budeme počítat ručně
data_ctyri_kvety

    sirka_listku delka_listku
3            0.2          1.3
82           1.0          3.7
117          1.8          5.5
101          2.5          6.0

TipKonkrétní kandidátní přímka

Jako konkrétní příklad použijeme přímku (B) s interceptem a = 1,1 a sklonem b = 2,2.

Residuum je naměřená délka minus odhadnutá délka v cm: svislý rozdíl mezi bodem a přímkou při stejné šířce. Úsečka znázorňuje jeho velikost, poloha bodu vůči přímce určuje znaménko. Pokud je kladné, model délku podhodnocuje; pokud záporné, délku nadhodnocuje.

Residua přitom očekáváme i u užitečného modelu. Jednotlivé organismy se v důsledku přirozené biologické variability liší, takže žádná jednoduchá přímka nepopíše každý květ přesně. Samotná existence residuí proto neznamená, že je model špatný.

V grafu níže jsou residua vidět přímo jako oranžové úsečky mezi očíslovanými šedými body a fialovou přímkou B.

5.2 Jaký rozdíl zůstane u jednoho květu?

Leží čtvrtý květ (vpravo nahoře) nad, nebo pod přímkou B? Z grafu odhadněte znaménko rozdílu; výpočtem jej ověříme.

Pro jeden konkrétní květ použijeme tento postup:

1. Vezmeme šířku vybraného korunního lístku.
2. Dosadíme ji do vztahu odhadnutá délka = intercept + sklon × šířka.
3. Vypočítáme odhadnutou délku tohoto korunního lístku podle zvolené přímky.

Odhadnutou délku vyjadřuje tento vztah mezi veličinami:

\[ \begin{aligned} \text{odhadnutá délka čtvrtého korunního lístku} &= \text{intercept} \\ &\quad + \text{sklon} \cdot \text{šířka čtvrtého korunního lístku}. \end{aligned} \]

Šířka čtvrtého květu 2,5 cm určuje vodorovnou polohu bodu 4 v grafu výše. Intercept 1,1 a sklon 2,2 popisují fialovou přímku. Po dosazení:

\[ \begin{aligned} \text{odhadnutá délka}_4 &= \text{intercept} + \text{sklon} \cdot \text{šířka}_4 \\ &= \text{1,1} + \text{2,2} \cdot \text{2,5} \\ &= \text{1,1} + \text{5,5} \\ &= \text{6,6}\ \text{cm} \end{aligned} \]

Výsledek 6,6 cm je v grafu výše výška fialové přímky při šířce 2,5 cm, tedy horní konec čtvrté oranžové úsečky. Naměřená délka 6 cm je naopak výška bodu 4. Residuum je rozdíl těchto dvou výšek:

\[ \begin{aligned} \text{residuum čtvrtého korunního lístku} &= \text{naměřená délka čtvrtého korunního lístku} \\ &\quad - \text{odhadnutá délka čtvrtého korunního lístku}. \end{aligned} \]

S hodnotami odečtenými z bodu a přímky:

\[ \begin{aligned} \text{residuum}_4 &= \text{naměřená délka}_4 - \text{odhadnutá délka}_4 \\ &= \text{6} - \text{6,6} \\ &= \text{−0,6}\ \text{cm} \end{aligned} \]

Hodnota -0,6 cm odpovídá čtvrté oranžové úsečce. Samotná úsečka má délku 0,6 cm; záporné znaménko residua říká, že bod 4 leží pod fialovou přímkou.

# Intercept (cm) a sklon (cm/cm) přímky B
a0 <- 1.1
b0 <- 2.2

# Dosadíme šířku čtvrtého květu do rovnice zvolené přímky
vec_odhadnuta_delka_4 <-
  a0 + b0 * data_ctyri_kvety$sirka_listku[4]

vec_odhadnuta_delka_4

[1] 6.6

5.3 Naměřená délka, odhad a residuum čtyř květů

Stejný výpočet nyní použijeme pro všechny čtyři květy. Získáme tak:

• odhadnuté délky jednotlivých květů
• rozdíly mezi naměřenou a odhadnutou délkou (residua)

vec_odhadnuta_delka <-
  a0 + b0 * data_ctyri_kvety$sirka_listku

# Residuum je naměřená hodnota minus odhad přímky
vec_residuum <-
  data_ctyri_kvety$delka_listku - vec_odhadnuta_delka

table_pred <-
  data.frame(
    kvet = seq_len(length.out = nrow(x = data_ctyri_kvety)),
    sirka = data_ctyri_kvety$sirka_listku,
    namerena_delka = data_ctyri_kvety$delka_listku,
    odhadnuta_delka = round(x = vec_odhadnuta_delka, digits = 2),
    residuum = round(x = vec_residuum, digits = 2)
  )

# Zobrazíme měření, odhad i residuum vedle sebe
table_pred

  kvet sirka namerena_delka odhadnuta_delka residuum
1    1   0.2            1.3            1.54    -0.24
2    2   1.0            3.7            3.30     0.40
3    3   1.8            5.5            5.06     0.44
4    4   2.5            6.0            6.60    -0.60

Čtení tabulky krok za krokem (všechny hodnoty jsou v cm):

• sloupec kvet: pořadí květu, stejné jako čísla bodů v grafu
• sloupec sirka: šířka korunního lístku konkrétního květu
• sloupec namerena_delka: skutečně naměřená délka korunního lístku
• sloupec odhadnuta_delka: délka odhadnutá pomocí přímky B
• sloupec residuum: rozdíl mezi naměřenou a odhadnutou délkou

Květy 2 a 3 mají kladné residuum (leží nad přímkou), květy 1 a 4 záporné (leží pod ní).

Stejný výpočet můžeme pro každý květ zapsat pomocí symbolů. Index \(i\) označuje pořadí květu. Naměřenou hodnotu označíme \(y_i\), šířku \(x_i\), odhadnutou hodnotu \(\hat{y}_i\), intercept \(a\), sklon \(b\) a residuum (anglicky residual) \(e_i\).

\[ \hat{y}_i = a + b \cdot x_i \]

\[ e_i = y_i - \hat{y}_i \]

V prvním vzorci je \(\hat{y}_i\) odhadnutá délka jednoho lístku, \(x_i\) jeho naměřená šířka, \(a\) intercept a \(b\) sklon. Ve druhém vzorci je \(y_i\) skutečně naměřená délka a \(e_i\) rozdíl mezi měřením a odhadem.

5.4 Jak porovnat vzdálenosti všech květů od přímky?

Nabízí se residua jednoduše sečíst. Pro čtyři květy a přímku B ale vyjde:

\[ \text{−0,24} + \text{0,40} + \text{0,44} + (\text{−0,60}) = \text{0}\ \text{cm}. \]

Kladná a záporná residua se vyruší, přestože přímka žádným ze čtyř bodů neprochází. Prostý součet proto nepozná, jak daleko body od přímky leží.

Součet čtverců residuí (anglicky sum of squared errors nebo residual sum of squares, zkratka SSE) tento problém řeší: každé residuum nejprve umocní na druhou. Říká, jak moc jsou body daleko od modelové přímky. Čím menší SSE, tím lépe přímka vystihuje data.

Proč čtverce? Protože:

• odstraní znaménko (kladná a záporná residua se neodečtou)
• více penalizují velké odchylky
• dávají jedno srovnatelné číslo pro každou kandidátní přímku

Výpočet slovy pro tuto přímku:

\[ \text{součet čtverců residuí} = \left(\text{residuum prvního květu}\right)^2 + \cdots + \left(\text{residuum posledního květu}\right)^2. \]

S residui z tabulky výše:

\[ \begin{aligned} \operatorname{SSE} &= (\text{−0,24})^2 + (\text{0,40})^2 \\ &\quad + (\text{0,44})^2 + (\text{−0,60})^2 \\ &= \text{0,0576} + \text{0,16} \\ &\quad + \text{0,1936} + \text{0,36} \\ &= \text{0,7712}\ \text{cm}^2. \end{aligned} \]

Obecně pro všechna pozorování:

\[ SSE = \sum_{i=1}^{n} e_i^2 \]

Zde \(e_i\) označuje residuum jednoho květu, \(e_i^2\) jeho čtverec, \(n\) počet květů a znak \(\sum\) znamená, že čtverce všech residuí sečteme.

Residua mají jednotku cm, jejich čtverce a SSE jednotku cm². Porovnávat SSE smíme pro stejné květy a stejné jednotky.

Stejný výpočet nyní provedeme pro každou kandidátní přímku a pro všech 150 květů. Výsledkem je jediné srovnatelné číslo pro každou přímku:

Porovnání kandidátních přímek ve všech datech

Přímka	Intercept (cm)	Sklon (cm/cm)	SSE (cm²)
A	0,5	1,8	231,08
B	1,1	2,2	33,98
C	1,7	2,6	214,32

Která z nabízených přímek má nejmenší SSE? Vysvětlete, proč její SSE nelze přímo porovnat s předchozím součtem pro čtyři květy.

DůležitéOdpověď: nejlepší přímka má nejmenší SSE

Z nabízených přímek je nejlepší B, protože má nejmenší SSE. Tabulka používá všech 150 květů, zatímco ruční příklad pouze čtyři; jednotlivé součty proto hodnotí jiné množství měření. Lineární model ale není omezen na tyto tři možnosti: hledá takový intercept a sklon, které SSE minimalizují mezi všemi možnými přímkami.

5.5 Naměřená délka se skládá z odhadu a residua

Naměřenou délku čtvrtého květu můžeme rozdělit na délku odhadnutou přímkou a residuum:

\[ \begin{aligned} \text{naměřená délka čtvrtého korunního lístku} &= \text{odhadnutá délka čtvrtého korunního lístku} \\ &\quad + \text{residuum čtvrtého korunního lístku}. \end{aligned} \]

S hodnotami čtvrtého květu:

\[ \begin{aligned} \text{naměřená délka}_4 &= \text{odhadnutá délka}_4 + \text{residuum}_4 \\ &= \text{6,6} + (\text{−0,60}) \\ &= \text{6}\ \text{cm} \end{aligned} \]

Pomocí již zavedených obecných symbolů:

\[ y_i = \hat{y}_i + e_i = a + b \cdot x_i + e_i \]

Interpretace v biologickém příběhu:

• \(y_i\) je naměřená délka jednoho korunního lístku a \(x_i\) jeho naměřená šířka
• \(\hat{y}_i\) je délka odhadnutá přímkou pro tuto šířku
• a (intercept) je odhad délky korunního lístku při nulové šířce; nula leží těsně za rozsahem pozorovaných dat
• b (sklon) říká, o kolik se v průměru liší odhadnutá délka dvou korunních lístků, jejichž šířky se liší o 1 cm
• \(e_i\) je nevysvětlená část variability (residuum); index \(i\) označuje pořadí květu

5.6 O kolik se liší odhadnutá délka?

Sklon zvolené přímky spočítáme ze dvou bodů, které na této přímce leží; nestačí dvě libovolná pozorování z dat. Použijeme šířky prvního a čtvrtého květu:

Dva body ležící na zvolené přímce

Bod na přímce	Šířka (cm)	Odhadnutá délka (cm)
při šířce květu 1	0,2	1,54
při šířce květu 4	2,5	6,60

Sklon přímky vyjadřuje vztah mezi dvěma změnami:

\[ \text{sklon přímky} = \frac{ \text{změna odhadnuté délky korunního lístku} }{ \text{změna šířky korunního lístku} }. \]

Se souřadnicemi dvou bodů na přímce:

\[ \text{sklon přímky} = \frac{\text{6,6} - \text{1,54}} {\text{2,5} - \text{0,2}} = \frac{\text{5,06}} {\text{2,3}} = \text{2,2} \]

Obecně:

\[ b = \frac{\Delta \hat{y}}{\Delta x} \]

Znak \(\Delta\) znamená změnu neboli rozdíl dvou hodnot. \(\Delta \hat{y}\) je rozdíl odhadnutých délek a \(\Delta x\) rozdíl šířek dvou bodů na přímce. Jednotka sklonu je cm délky na cm šířky (cm/cm); obě veličiny jsou délky, ale jejich biologické role zůstávají odlišné.

6 Přímka pro délku korunního lístku v R

Funkce lm() fituje lineární model. Její zápis si přečteme jako obecnou gramatiku:

# Obecný zápis pojmenuje role proměnných a dat
lm(
  formula = odezva ~ prediktor,
  data = data
)

V tomto zápisu jsou odezva, prediktor a data zástupné názvy. Znak ~ můžeme číst jako „odezvu vysvětlujeme pomocí prediktoru“.

Pro naši biologickou otázku dosadíme:

• odezva: delka_listku
• prediktor: sirka_listku
• data: data_kosatce

Konkrétní volání odděluje dvě operace:

• lm() odhadne parametry modelu z dat
• coef(object = mod_listky) vrátí dva parametry, které už znáte z kandidátních přímek

# Odhadneme vztah mezi šířkou a délkou korunního lístku
mod_listky <-
  lm(
    formula = delka_listku ~ sirka_listku,
    data = data_kosatce
  )

# Zobrazíme odhadnutý intercept a sklon
coef(object = mod_listky)

 (Intercept) sirka_listku 
    1.083558     2.229940 

Co vidíme za parametry (koeficienty) modelu:

• Intercept odpovídá parametru a
• odhad u sirka_listku odpovídá sklonu b

Přečtěte odhad u sirka_listku v biologických jednotkách. Jak se podle modelu liší odhadnuté délky dvou květů, jejichž šířky se liší o 1 cm? Platí tento rozdíl přesně pro každou dvojici květů?

DůležitéOdpověď na biologickou otázku

Odhadnutý sklon je 2,23 cm/cm. Když porovnáme dva květy, jejichž korunní lístky se liší šířkou o 1 cm, model pro květ se širším lístkem odhaduje v průměru o 2,23 cm delší korunní lístek.

Pro každou konkrétní dvojici květů tento rozdíl přesně neplatí. Jednotlivé květy mají residua, takže rozdíl jejich naměřených délek může být větší i menší. Sklon popisuje průměrný vztah, ne pravidlo pro každý květ.

Odhadnutý intercept je 1,08 cm. Popisuje odhadnutou délku při nulové šířce, která leží těsně mimo rozsah pozorovaných dat. V tomto biologickém příběhu je proto důležitější sklon.

Model tak odpovídá na směr a velikost vztahu. Ani dobře vypadající lineární vztah ale automaticky neznamená kauzální vztah: z modelu zatím čteme pouze asociaci mezi proměnnými.

6.1 Přímka odhadnutá ze všech květů

Odhadnuté délky pro všechny květy vypočítáme pomocí predict(). Je to stejný výpočet jako ruční dosazení u čtyř květů, jen pro celý dataset:

vec_odhadnuta_delka <-
  predict(
    object = mod_listky,
    newdata = data_kosatce
  )

# Zkontrolujeme prvních šest odhadnutých délek
head(x = vec_odhadnuta_delka)

       1        2        3        4        5        6 
1.529546 1.529546 1.529546 1.529546 1.529546 1.975534 

Graf naměřených hodnot s přímkou odhadnutou modelem ukazuje, že přímka opravdu odpovídá číslům z modelu:

Základní graf v R

# Naměřené dvojice zůstávají body ve svých původních jednotkách
plot(
  x = data_kosatce$sirka_listku,
  y = data_kosatce$delka_listku,
  xlab = "Šířka korunního lístku (cm)",
  ylab = "Délka korunního lístku (cm)",
  pch = 16,
  col = "#8A8A8A",
  cex.axis = 1.4,
  cex.lab = 1.4
)

# Přímku určuje intercept a sklon odhadnutého modelu
abline(
  a = coef(object = mod_listky)[["(Intercept)"]],
  b = coef(object = mod_listky)[["sirka_listku"]],
  col = "#5D2890",
  lwd = 2
)

Graf pomocí {ggplot2}

data_model <-
  data.frame(
    # Původní naměřené hodnoty
    sirka_listku = data_kosatce$sirka_listku,
    delka_listku = data_kosatce$delka_listku,
    # Délky odhadnuté stejným modelem
    odhadnuta_delka = predict(
      object = mod_listky,
      newdata = data_kosatce
    )
  )

# Body jsou naměřené květy; čára spojuje odhady modelu
ggplot2::ggplot(
  data = data_model,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.75,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_line(
    ggplot2::aes(y = odhadnuta_delka),
    color = "#5D2890",
    linewidth = 1
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::theme_bw(base_size = 18)

TipDoplňující: kratší alternativa pomocí geom_smooth()

Balíček {ggplot2} umí lineární model odhadnout a nakreslit přímo pomocí geom_smooth(method = "lm"). Tuto zkratku používejte hlavně pro první rychlou vizualizaci vztahu. Vrstva sice model interně odhadne, ale nevytvoří snadno dostupný objekt modelu, ze kterého bychom mohli získat koeficienty, residua nebo další výsledky. Pro vlastní analýzu proto model fitujeme zvlášť pomocí lm() a uložíme jej do objektu. Argument se = FALSE zde vypíná zobrazení pásu nejistoty, ke kterému se vrátíme v další lekci.

ggplot2::ggplot(
  # Použijeme původní měření šířky a délky
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.75,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_smooth(
    # Funkce odhadne lineární model a rovnou nakreslí jeho přímku
    method = "lm",
    se = FALSE,
    color = "#5D2890"
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::theme_bw(base_size = 18)

7 Co zůstává mezi měřením a odhadem?

Residua jsou uložena přímo v objektu modelu. Získáme je funkcí resid() a uložíme do vektoru, se kterým můžeme dále pracovat:

vec_rezidua <-
  # resid() vrátí residuum pro každý květ použitý při fitování
  resid(object = mod_listky)

# Zkontrolujeme prvních šest residuí
head(x = vec_rezidua)

          1           2           3           4           5           6 
-0.12954613 -0.12954613 -0.22954613 -0.02954613 -0.12954613 -0.27553423 

Kladné residuum znamená, že naměřená délka leží nad odhadem modelu. Záporné residuum znamená, že naměřená délka leží pod odhadem modelu.

7.1 Jak jsou residua rozdělená?

Histogram už znáte z předchozích lekcí. Která znaménka jsou častá a jak daleko residua sahají od nuly? Tentokrát na vodorovné ose nejsou původní délky lístků, ale rozdíly mezi naměřenou a odhadnutou délkou:

hist(
  # Histogram ukazuje rozdělení residuí
  x = vec_rezidua,
  breaks = 12,
  main = "Histogram residuí",
  xlab = "Residuum (cm)",
  ylab = "Četnost",
  col = "#D9D9D9",
  border = "white",
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Svislá čára označí nulové residuum
  v = 0,
  lty = 2,
  col = "#5D2890"
)

Přerušovaná čára označuje nulové residuum, tedy přesnou shodu naměřené a odhadnuté délky. Histogram ukazuje, která residua jsou častá a zda některé květy leží od odhadu neobvykle daleko. Neříká ale, u kterých odhadnutých délek tato residua vznikají.

7.2 Mění se residua podél odhadnuté délky?

Proto ke každému květu spojíme odhadnutou délku z fitted() s jeho residuem z resid():

vec_odhadnuta_delka <-
  # fitted() vrátí odhadnutou délku každého květu
  fitted(object = mod_listky)

data_rezidua <-
  data.frame(
    odhadnuta_delka = vec_odhadnuta_delka,
    reziduum = vec_rezidua
  )

# Zkontrolujeme, že každý odhad má odpovídající residuum
head(x = data_rezidua)

  odhadnuta_delka    reziduum
1        1.529546 -0.12954613
2        1.529546 -0.12954613
3        1.529546 -0.22954613
4        1.529546 -0.02954613
5        1.529546 -0.12954613
6        1.975534 -0.27553423

Stejný graf ukazují dvě záložky: v základním R a pomocí balíčku {ggplot2}.

Základní graf v R

plot(
  # Každý bod spojuje odhad modelu s residuem stejného květu
  x = data_rezidua$odhadnuta_delka,
  y = data_rezidua$reziduum,
  xlab = "Odhadnutá délka (cm)",
  ylab = "Residuum (cm)",
  main = "Odhadnuté délky vs. residua",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # V dobře padnoucím modelu hledáme body rozptýlené kolem nuly
  h = 0,
  lty = 2,
  col = "#5D2890"
)

Graf pomocí {ggplot2}

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_rezidua,
  mapping = ggplot2::aes(x = odhadnuta_delka, y = reziduum)
) +
  ggplot2::geom_hline(
    # Referenční čára ukazuje nulové residuum
    yintercept = 0,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Odhadnutá délka (cm)",
    y = "Residuum (cm)",
    title = "Odhadnuté délky vs. residua"
  ) +
  ggplot2::theme_bw(base_size = 18)

Jak tento graf číst:

• každý bod představuje jeden květ
• vodorovná poloha ukazuje délku odhadnutou modelem
• svislá poloha ukazuje rozdíl mezi naměřenou a odhadnutou délkou
• bod na nulové čáře má naměřenou délku přesně rovnou odhadu
• bod nad nulou model podhodnotil, bod pod nulou model nadhodnotil

Pro první kontrolu stačí jednoduchá otázka: tvoří body kolem nuly neuspořádaný mrak, nebo v nich zůstává viditelný vzor? Viditelný vzor znamená, že model v datech něco systematicky nezachytil.

Najděte jednu oblast grafu, kde residua leží převážně nad nebo pod nulou. Co to říká o odhadovaných délkách v této oblasti? Směr odchylky pojmenujte bez návrhu nového modelu; k biologickému vysvětlení se vrátíme s informací o druhu.

8 Co musí být rozumné, aby přímka odpovídala na otázku?

Lineární model je zjednodušený popis dat. Nemusí být dokonale pravdivý, ale jeho zjednodušení nesmí být v přímém rozporu s otázkou a daty. U jednoduchého modelu sledujeme čtyři odlišné oblasti:

1. Tvar vztahu: průměrná délka by se v pozorovaném rozsahu šířek měla měnit přibližně po přímce. Zřetelný oblouk v residuích by upozornil, že lineární tvar nestačí.
2. Rozptýlení kolem přímky: svislý rozptyl residuí by se neměl systematicky zvětšovat nebo zmenšovat. Trychtýřovitý tvar by znamenal, že model zachází se všemi částmi vztahu jako se stejně proměnlivými, přestože data tak nevypadají.
3. Tvar residuí: aby pozdější odhady nejistoty byly přiměřené, nesmí tvar residuí ukazovat silný rozpor s použitým modelem. Několik drobných odchylek od ideálního tvaru je v reálných datech běžných.
4. Nezávislost pozorování: zda jsou pozorování nezávislá, vyplývá především ze způsobu získání dat. Graf residuí sám nepozná, že několik měření pochází ze stejného jedince, lokality nebo časové řady.

8.1 Čtyři otázky, čtyři různé stopy v grafu

Následující schémata nejsou další čtyři datasety ani automatický test. Jsou to čtyři záměrně zjednodušené obrazy, které pomáhají pojmenovat, na jaký druh problému se díváme.

Schématické stopy v grafech residuí. Vodorovná čárkovaná čára označuje nulové residuum; spodní pravý panel jako jediný používá na vodorovné ose pořadí měření.

Vlevo nahoře nevidíme nápadný rozpor, ale ani tento obraz správnost modelu nedokazuje. Oblouk upozorňuje na tvar vztahu, trychtýř na měnící se rozptýlení. Dlouhé běhy nad a pod nulou se ukážou teprve tehdy, když zachováme pořadí měření; nezávislost proto posuzujeme společně s informací o tom, jak data vznikla.

8.2 Odezva a residua mohou mít jiný tvar

U kosatců můžeme porovnat dvě skutečná rozdělení. Levý histogram popisuje naměřené délky korunních lístků. Pravý histogram popisuje odchylky těchto délek od hodnot odhadnutých přímkou.

Naměřená odezva a residua stejného modelu odpovídají na jiné otázky, a proto jejich rozdělení nemusí mít stejný tvar. Vodorovné osy panelů mají odlišný rozsah.

Histogram odezvy se ptá, jaké délky jsme naměřili. Histogram residuí se ptá, jak se skutečné délky liší od přímky při dané šířce. Ani jeden histogram sám neposoudí linearitu, měnící se rozptýlení nebo nezávislost; k tomu potřebujeme předchozí grafy a znalost uspořádání studie.

Diagnostika proto není zkouška, kterou model jednou provždy „udělá“ nebo „neudělá“. Spojuje grafy, znalost uspořádání studie a otázku, zda případný problém může změnit biologický závěr. V této lekci problémy pouze rozpoznáváme; jejich řešení přijdou postupně v dalších lekcích.

TipDoplňující: předpoklady se týkají residuí kolem modelu

Lineární model nevyžaduje, aby samotná šířka nebo délka korunních lístků měla dokonale normální rozdělení. Důležité je, jak se chovají odchylky délky od hodnot očekávaných při dané šířce. Proto diagnostikujeme residua kolem přímky, nikoli pouze histogram původní odezvy.

Ani residua nemusí vytvořit dokonalou zvonovitou křivku. U malého souboru lze tvar posoudit jen hrubě; u velkého souboru zase i nepatrná odchylka může vypadat „statisticky zjistitelně“. Smyslem kontroly je hledat odchylku natolik výraznou, že zpochybňuje použitý model nebo jeho nejistotu, nikoli získat mechanické razítko správnosti.

TipDoplňující: Q–Q graf residuí

Tvar histogramu residuí závisí také na tom, kolik sloupců zvolíme. Podrobnější pohled na tvar residuí nabízí Q–Q graf (z anglického quantile–quantile). Využívá kvantily – obecný název pro percentily a kvartily z první lekce, tedy hodnoty, pod nimiž leží určitý podíl seřazených dat. Residua seřadíme od nejmenšího po největší a každé porovnáme s hodnotou, kterou bychom na stejném místě pořadí čekali u souměrného zvonovitého (normálního) rozdělení.

Základní graf v R

qqnorm(
  # Seřazená residua porovnáme s kvantily normálního rozdělení
  y = vec_rezidua,
  main = "Q–Q graf residuí",
  xlab = "Teoretické kvantily normálního rozdělení",
  ylab = "Residuum (cm)",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

qqline(
  # Referenční přímka prochází prvním a třetím kvartilem residuí
  y = vec_rezidua,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)

Graf pomocí {ggplot2}

ggplot2::ggplot(
  # Stejná residua jako v předchozím grafu v základním R
  data = data_rezidua,
  mapping = ggplot2::aes(sample = reziduum)
) +
  ggplot2::stat_qq_line(
    # Referenční přímka prochází prvním a třetím kvartilem residuí
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::stat_qq(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Teoretické kvantily normálního rozdělení",
    y = "Residuum (cm)",
    title = "Q–Q graf residuí"
  ) +
  ggplot2::theme_bw(base_size = 18)

Jak tento graf číst:

• každý bod je residuum jednoho květu; svislá osa ukazuje jeho skutečnou hodnotu v cm
• vodorovná osa ukazuje, kde by residuum se stejným pořadím leželo u normálního rozdělení, měřeno ve směrodatných odchylkách od středu
• přerušovaná přímka ukazuje, kde by body ležely při dokonale zvonovitém tvaru
• body blízko přímky: tvar residuí se od zvonovitého rozdělení výrazně neliší
• oba konce se odklánějí opačnými směry (vlevo dolů, vpravo nahoru): velkých kladných i záporných residuí je víc, než by odpovídalo zvonovité křivce
• body tvoří oblouk: residua jsou nesouměrná, na jedné straně nuly sahají dál než na druhé

U kosatců leží prostřední residua téměř přesně na přímce. Na obou koncích se body mírně odklánějí opačnými směry: několik květů má o něco větší kladná i záporná residua, než by odpovídalo zvonovité křivce. Odchylka se týká hlavně nejkrajnějších bodů a je mírná; sama o sobě proto použití lineárního modelu nezpochybňuje.

Q–Q graf není test, který model „schválí“. I residua vzniklá přesně ze zvonovitého rozdělení se na koncích grafu trochu vlní, zvlášť u malého počtu pozorování. Hledáme proto výraznou a souvislou odchylku, ne každý bod mimo přímku.

TipDoplňující: odlehlý, krajní a vlivný bod

Odlehlé pozorování má neobvykle velkou odchylku od toho, co jsme čekali. Krajní bod má neobvyklou hodnotu prediktoru, například mimo většinu pozorovaných šířek. Vlivný bod je takový, jehož přítomnost výrazně mění odhadnutý intercept nebo sklon.

Jeden květ může splňovat jednu, dvě nebo všechny tři vlastnosti. Žádná z nich sama neznamená chybu ani povolení řádek smazat. Nejdříve ověříme měření a biologický kontext, potom zkoumáme citlivost závěru. Jak „krajnost“ a vliv jednotlivých květů změřit, ukazují dva následující doplňující rámečky.

TipDoplňující: pákový efekt

Pákový efekt (anglicky leverage) měří, jak neobvyklou hodnotu prediktoru má jedno pozorování. Přímku si můžeme představit jako páku: květ, jehož šířka leží daleko od průměrné šířky, ji může naklonit snáz než květ uprostřed rozsahu. V R pákový efekt získáme funkcí hatvalues():

vec_pakovy_efekt <-
  # hatvalues() vrátí pákový efekt pro každý květ
  hatvalues(model = mod_listky)

# Spojíme měření, residua a pákový efekt jednotlivých květů
data_vliv_kvetu <-
  data.frame(
    sirka_listku = data_kosatce$sirka_listku,
    delka_listku = data_kosatce$delka_listku,
    reziduum = vec_rezidua,
    pakovy_efekt = vec_pakovy_efekt
  )

# Zkontrolujeme nejmenší a největší pákový efekt
range(vec_pakovy_efekt)

[1] 0.006666672 0.026208483

Pákový efekt nemá jednotky. S jediným prediktorem závisí pouze na šířce lístku: čím dál je šířka od průměrné šířky, tím je hodnota větší, a to na obě strany. Nejmenší hodnotu proto mají květy se šířkou blízko průměru. Naměřená délka ani residuum na pákový efekt vliv nemají. V následujícím grafu proto vyneseme obě vlastnosti z předchozího rámečku proti sobě: krajnost na vodorovnou osu a odlehlost na svislou.

Základní graf v R

plot(
  # Vodorovně krajnost šířky, svisle odchylka od přímky
  x = data_vliv_kvetu$pakovy_efekt,
  y = data_vliv_kvetu$reziduum,
  xlab = "Pákový efekt",
  ylab = "Residuum (cm)",
  main = "Pákový efekt vs. residua",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Nulové residuum: naměřená délka je rovna odhadu
  h = 0,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)

Graf pomocí {ggplot2}

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_vliv_kvetu,
  mapping = ggplot2::aes(x = pakovy_efekt, y = reziduum)
) +
  ggplot2::geom_hline(
    # Nulové residuum: naměřená délka je rovna odhadu
    yintercept = 0,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Pákový efekt",
    y = "Residuum (cm)",
    title = "Pákový efekt vs. residua"
  ) +
  ggplot2::theme_bw(base_size = 18)

Jak tento graf číst:

• vpravo leží květy s nejvíce neobvyklou šířkou lístku
• daleko od nulové čáry leží květy, jejichž délka se nejvíce liší od odhadu
• body tvoří svislé sloupce: mnoho květů má stejnou šířku a stejný pákový efekt mají i květy stejně daleko od průměrné šířky na opačných stranách – velmi úzké i velmi široké lístky

Úplně vpravo leží květy s nejširšími lístky (2,5 cm). Všechny leží pod nulovou čarou: model jejich délku nadhodnocuje. Květ, který je zároveň vpravo a daleko od nuly, může přímku naklonit nejvíce. Právě tuto kombinaci shrnuje Cookova vzdálenost.

TipDoplňující: Cookova vzdálenost

Cookova vzdálenost (anglicky Cook’s distance) spojuje obě vlastnosti z předchozího rámečku do jednoho čísla. Pro každý květ vyjadřuje, jak moc by se změnily odhadnuté délky všech květů, kdybychom právě tento květ vynechali a přímku odhadli znovu. Velkou hodnotu mívá hlavně květ, který má velké residuum a zároveň velký pákový efekt.

vec_cookova_vzdalenost <-
  # cooks.distance() vrátí Cookovu vzdálenost pro každý květ
  cooks.distance(model = mod_listky)

# Doplníme Cookovu vzdálenost a pořadí květu v datech
data_vliv_kvetu$cookova_vzdalenost <- vec_cookova_vzdalenost
data_vliv_kvetu$poradi_kvetu <-
  seq_len(length.out = nrow(x = data_vliv_kvetu))

# Orientační hodnota 4/n, kde n je počet květů v modelu
hranice_cook <- 4 / nrow(x = data_vliv_kvetu)
hranice_cook

[1] 0.02666667

Základní graf v R

plot(
  # Každá svislá čárka představuje jeden květ
  x = data_vliv_kvetu$poradi_kvetu,
  y = data_vliv_kvetu$cookova_vzdalenost,
  type = "h",
  lwd = 2,
  col = "#8A8A8A",
  xlab = "Pořadí květu v datech",
  ylab = "Cookova vzdálenost",
  main = "Cookova vzdálenost květů",
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Orientační čára, nikoli hranice pro mazání dat
  h = hranice_cook,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)

Graf pomocí {ggplot2}

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_vliv_kvetu,
  mapping = ggplot2::aes(x = poradi_kvetu, y = cookova_vzdalenost)
) +
  ggplot2::geom_segment(
    # Každá svislá čárka vede od nuly k hodnotě jednoho květu
    mapping = ggplot2::aes(xend = poradi_kvetu, yend = 0),
    linewidth = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_hline(
    # Orientační čára, nikoli hranice pro mazání dat
    yintercept = hranice_cook,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::labs(
    x = "Pořadí květu v datech",
    y = "Cookova vzdálenost",
    title = "Cookova vzdálenost květů"
  ) +
  ggplot2::theme_bw(base_size = 18)

Nad orientační čarou leží 11 z 150 květů. Hodnota 4/n není test ani pravidlo pro mazání dat; pouze upozorňuje, které květy si zaslouží bližší pohled. Seřadíme je od největší Cookovy vzdálenosti:

data_vliv_serazene <-
  data_vliv_kvetu[
    # Seřadíme květy od největší Cookovy vzdálenosti
    order(
      data_vliv_kvetu$cookova_vzdalenost,
      decreasing = TRUE
    ),
  ]

# Vybereme sloupce, které chceme u květů porovnat
vec_sloupce_tabulky <-
  c(
    "sirka_listku",
    "delka_listku",
    "reziduum",
    "pakovy_efekt",
    "cookova_vzdalenost"
  )

# Zobrazíme pět nejvlivnějších květů se zaokrouhlenými hodnotami;
# čísla vlevo jsou pořadí květů v datech
round(
  x = head(x = data_vliv_serazene[, vec_sloupce_tabulky], n = 5),
  digits = 3
)

    sirka_listku delka_listku reziduum pakovy_efekt cookova_vzdalenost
115          2.4          5.1   -1.335        0.023              0.095
142          2.3          5.1   -1.112        0.021              0.058
145          2.5          5.7   -0.958        0.026              0.056
146          2.3          5.2   -1.012        0.021              0.048
123          2.0          6.7    1.157        0.014              0.042

Nejvlivnější květ má široký lístek (2,4 cm) a jeho naměřená délka leží 1,34 cm pod odhadem přímky. Ke kterým druhům patří květy nad orientační čarou?

# Spočítáme druhy květů s Cookovou vzdáleností nad orientační čarou
table(data_kosatce$druh[vec_cookova_vzdalenost > hranice_cook])

virginica 
       11 

Všechny tyto květy patří druhu virginica, který má nejširší lístky, a tedy i velký pákový efekt. Model 8 z nich nadhodnocuje a 3 naopak výrazně podhodnocuje. Proč jedna společná přímka pro tři druhy nestačí, ukazuje kapitola Které druhy model nadhodnocuje a podhodnocuje?.

Jak moc nejvlivnější květ přímkou skutečně hýbe, zjistíme přímo: odhadneme stejný model bez něj a porovnáme koeficienty.

# Najdeme řádek květu s největší Cookovou vzdáleností
radek_nejvlivnejsi <- which.max(x = vec_cookova_vzdalenost)

mod_bez_kvetu <-
  lm(
    # Stejný model odhadneme bez jediného květu
    formula = delka_listku ~ sirka_listku,
    data = data_kosatce[-radek_nejvlivnejsi, ]
  )

# Porovnáme intercept a sklon obou modelů
rbind(
  vsechny_kvety = coef(object = mod_listky),
  bez_nejvlivnejsiho = coef(object = mod_bez_kvetu)
)

                   (Intercept) sirka_listku
vsechny_kvety         1.083558     2.229940
bez_nejvlivnejsiho    1.069930     2.248904

Sklon se změní z 2,23 na 2,25 cm délky na 1 cm šířky. Ani nejvlivnější květ tedy přímku podstatně neposune a závěr o vztahu šířky a délky se nemění. Ani vysoká Cookova vzdálenost ale sama o sobě není důvod řádek smazat. Je to důvod ověřit měření, zvážit biologický kontext a porovnat závěr s květem a bez něj – tak jako zde.

9 Které druhy model nadhodnocuje a podhodnocuje?

Ve všech dosavadních grafech vypadaly květy stejně a druh kosatce v modelu není. Barvu teď použijeme pouze jako pomůcku, abychom zjistili, zda se jednotlivé druhy v datech chovají podobně.

9.1 Druhy v původních měřeních

Stejné barvy použijeme pro každý druh v obou následujících grafech:

Porovnejte polohu tří barevných skupin. Shrnuje společná přímka pouze změny uvnitř jednoho druhu, nebo také rozdíly mezi druhy?

Barvy ukazují, že tři druhy netvoří jeden společný bodový mrak. Jedna přímka vedená přes všechny květy proto shrnuje nejen vztah šířky a délky, ale také rozdíly mezi druhy.

9.2 Druhy v grafu odhadnutých délek a residuí

Ke každému residuu doplníme druh stejného květu:

Které barevné skupiny leží častěji nad nulou a které pod ní? Pro každý druh rozhodněte, zda model jeho délky spíše podhodnocuje, nebo nadhodnocuje. Použijte stejnou konvenci naměřená minus odhadnutá jako u jednoho květu.

DůležitéOdpověď: model jednotlivé druhy odhaduje různě

U versicolor leží nad nulou 84 % residuí, takže model jejich délky spíše podhodnocuje. U setosa leží pod nulou 78 % residuí a model jejich délky spíše nadhodnocuje. U virginica leží pod nulou 64 % residuí; nadhodnocení se týká hlavně květů s nejširšími lístky vpravo v grafu.

Po obarvení vidíme, že residua jednotlivých druhů neleží kolem nuly úplně stejně. Model některé druhy častěji nadhodnocuje a jiné častěji podhodnocuje. Informace o druhu tedy zůstala v části dat, kterou jednoduchý model nevysvětlil.

V této lekci stačí toto omezení rozpoznat a popsat. Druh jako další prediktor přidáme do modelu až v některé z navazujících lekcí.

10 Co přímka popisuje a co zůstává otevřené?

Jednoduchý lineární model je výborný začátek, ale ne univerzální řešení. V této lekci jsme odhalili jeden důvod k opatrnosti: v residuích zůstaly rozdíly mezi druhy kosatců.

Z grafů již umíme rozpoznat varovné vzory:

• zakřivený průběh může ukazovat na nelinearitu
• postupné rozšiřování nebo zužování bodového mraku souvisí s pojmem heteroskedasticita
• jednotlivé body daleko od ostatních mohou být odlehlá pozorování
• barevně oddělené skupiny mohou upozornit na důležitý prediktor, který v modelu chybí

V této lekci stačí vzor graficky rozpoznat a popsat. Formální testy ani nápravu modelu zatím nevyžadujeme. Později podle biologické otázky zvážíme další prediktor, transformaci, jiný tvar vztahu nebo jiný model.

Hlavní princip zůstává stejný: nejprve omezení modelu transparentně pojmenujeme a teprve potom volíme vhodné rozšíření.

11 Co bude následovat dál

V následující lekci ponecháme stejný modelový rámec a přidáme první inferenční vrstvu: výběrovou variabilitu a intervaly spolehlivosti jako vyjádření nejistoty odhadu. Testování hypotéz a p-hodnota přijdou až v pozdější lekci.

12 Shrnutí

• Biologickou otázku převádíme na vztah mezi odezvou a prediktorem; u kosatců je délka korunního lístku odezva a šířka prediktor.
• Funkce lm() najde přímku s nejmenším součtem čtverců residuí. Její sklon udává, o kolik se v průměru liší odezva dvou pozorování, jejichž prediktor se liší o jednu jednotku.
• Residuum je rozdíl mezi naměřenou a odhadnutou hodnotou; residua nejsou „odpad“, ale informace o tom, co model nevysvětluje.
• V grafu residuí hledáme neuspořádaný mrak kolem nuly; zbývající vzor, například oddělené skupiny druhů, upozorňuje na omezení modelu.

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

Někdo vám ukáže jen odhadnutý sklon 2,23 cm/cm pro šířku a délku korunních lístků kosatců, bez grafu, a napíše titulek „O centimetr širší lístek, o 2,23 cm delší lístek“. Jaké tři otázky mu položíte, než titulku uvěříte?

DůležitéMožné otázky

Jak vypadá bodový graf a graf residuí — tvoří residua neuspořádaný mrak kolem nuly, nebo v nich zůstává vzor? Jsou v datech skupiny, například druhy, které vztah vytvářejí nebo mění? Kolik pozorování data obsahují a jak byla měřena? Jak přesný je odhad sklonu? A tvrdí titulek jen průměrný rozdíl mezi květy, nebo i to, že širší lístek délku způsobuje?

Slovníček pojmů

Pojmy použité v těchto materiálech a jejich definice. Úplný průběžně aktualizovaný slovníček je dostupný také online.

term	definition
bodovy-graf	Graf vztahu dvou numerických proměnných, v němž každý bod představuje jedno pozorování a jeho dvě hodnoty určují polohu na osách.
cookova-vzdalenost	Míra vlivu jednoho pozorování na model: vyjadřuje, jak moc by se změnily odhadnuté hodnoty, kdybychom toto pozorování vynechali.
diagnostika-modelu	Soubor grafických a numerických kontrol ověřujících, zda jsou splněny předpoklady modelu: linearita, normalita reziduí, homoskedasticita, absence vlivných bodů.
efekt	Kvantitativní popis toho, jak se podle modelu mění odezva s prediktorem; zahrnuje směr, velikost a nejistotu odhadu, ale sám o sobě nedokládá příčinu.
heteroskedasticita	Jev, kdy rozptyl reziduí není konstantní podél hodnot prediktoru - typicky se rozptyl zvyšuje s rostoucí predikovanou hodnotou.
histogram	Graf zobrazující rozdělení hodnot numerické proměnné: hodnoty jsou seřazeny do tříd, výška sloupce udává četnost v dané třídě.
intercept	Parametr modelu: hodnota odezvy predikovaná modelem, když jsou všechny prediktory rovno nule; geometricky jde o úsek na ose y (kde regresní přímka protíná osu y).
interval-spolehlivosti	Rozsah hodnot, v němž se s danou spolehlivostí (typicky 95 %) nachází odhadovaný parametr.
kategoricka-promenna	Proměnná nabývající hodnot z omezeného počtu kategorií; může být nominální (bez pořadí) nebo ordinální (s pořadím).
korelace	Bezrozměrné číslo popisující směr a těsnost vztahu dvou proměnných; Pearsonova korelace shrnuje lineární vztah a Spearmanova monotónní vztah podle pořadí.
kvantil	Hodnota, pod níž leží určený podíl seřazených dat; medián je 50% kvantil.
linearni-model	Statistický model předpokládající, že odezva je lineární funkcí prediktorů plus chybový člen s normálním rozdělením (y = a + b·x + ε).
median	Hodnota rozdělující seřazená data na dvě stejně početné poloviny; 50. percentil.
nejistota	Rozsah, v němž nevíme přesnou hodnotu parametru kvůli omezenému a variabilnímu vzorku; kvantifikuje se pomocí intervalů spolehlivosti nebo směrodatných chyb.
nelinearita	Situace, kdy vztah mezi prediktorem a odezvou není dobře popsatelný přímkou - má zakřivení nebo jinou nelineární formu.
numericka-promenna	Proměnná, jejíž hodnoty jsou čísla se smysluplnou velikostí nebo rozdílem; může být diskrétní, například počet, nebo spojitá, například délka.
odezva	Proměnná, jejíž variabilitu se snažíme vysvětlit nebo předpovědět; v notaci modelu `y ~ x` je odezva `y`.
odhadnuta-hodnota	Hodnota odezvy, kterou fitovaný model přiřadí pozorování na základě jeho hodnot prediktorů; značí se často ŷ.
odlehla-pozorovani	Hodnota, která se výrazně odlišuje od ostatních pozorování; může být způsobena chybou měření nebo jít o skutečně extrémní případ.
p-hodnota	Pravděpodobnost, že bychom při platnosti nulové hypotézy pozorovali data alespoň tak extrémní jako naše skutečná data.
pakovy-efekt	Míra toho, jak neobvyklé hodnoty prediktorů má jedno pozorování ve srovnání s ostatními; nezávisí na naměřené odezvě.
parametry-modelu	Číselné hodnoty odhadnuté z dat definující tvar modelu; v lineárním modelu jsou to intercept a sklon(y).
pozorovani	Jeden zaznamenaný případ v datech, který obsahuje hodnoty jedné nebo více proměnných.
prediktor	Proměnná, pomocí níž vysvětlujeme variabilitu odezvy; v notaci `y ~ x` je prediktor `x`.
promenna	Vlastnost zaznamenaná u pozorování, která může mezi pozorováními nabývat různých hodnot, například druh, délka nebo počet.
qq-graf	Graf, který porovnává seřazené hodnoty (typicky residua) s kvantily očekávanými u zvoleného rozdělení, nejčastěji normálního.
residuum	Rozdíl mezi pozorovanou hodnotou odezvy a hodnotou predikovanou modelem: eᵢ = yᵢ − ŷᵢ.
rozdeleni	Popis toho, jaké hodnoty proměnná nabývá a jak často se vyskytují; může jít o pozorované rozdělení dat nebo o teoretické pravděpodobnostní rozdělení.
rozptyl	Průměr čtverců odchylek od průměru; míra variability dat. Odmocnina rozptylu je směrodatná odchylka.
sklon	Parametr lineárního modelu udávající, o kolik se změní odezva, pokud prediktor vzroste o jednu jednotku (při zachování ostatních prediktorů).
smerodatna-odchylka	Odmocnina rozptylu; vyjadřuje průměrnou vzdálenost hodnot od průměru ve stejných jednotkách jako původní data.
soucet-ctvercu-residui	Součet druhých mocnin všech residuí; vyjadřuje celkovou čtvercovou odchylku pozorovaných hodnot od hodnot odhadnutých modelem.
transformace	Matematická operace použitá na hodnoty proměnné s cílem splnit předpoklady modelu nebo zlepšit linearitu vztahu.
variabilita	Obecný termín pro rozsah, v němž se hodnoty proměnné liší od sebe navzájem; kvantifikuje se pomocí rozptylu, SD, IQR nebo rozsahu.
vyberova-variabilita	Proměnlivost odhadu mezi možnými srovnatelnými výběry získanými stejným postupem.