Jak poznáme, že je model dobrý?

Logo kurzu Biostatistika BIOSTATISTIKA

L09

Porovnání kandidátních modelů: R², adjustované R² a AIC

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

Liška a veverka u terénního stolku s pásmem pod třemi různými jehličnany; opřené dřevěné destičky ukazují tři uspořádání fialových tyček a sova je z větve zamyšleně posuzuje. Ilustrace vytvořená pomocí AI.

Ilustrace vytvořená pomocí AI.

Co vlastně na stromě měříme?

Liška na lesní měřicí stanici omotává pásmo kolem kmene ve výšce prsou a veverka zaměřuje vrchol stromu; vedle sebe stojí tři různé jehličnany.

Ilustrace vytvořená pomocí AI.

U každého jehličnanu známe

  • druh stromu — tři severoamerické jehličnany;
  • průměr kmene ve výšce 1,3 m (cm);
  • výšku stromu (m).

Mají stejně silné stromy různých druhů stejnou výšku?

Co si pamatujete z minulé lekce?

Odpovězte samostatně na tři krátké otázky. Na telefonu neuvidíte odpovědi ostatních; společné výsledky otevřeme až po hlasování.

QR kód otevírající anonymní dotazník PollsLive se třemi otázkami.

Otevřít dotazník v novém okně

Jak odpovídala skupina?

Otevřít výsledky mimo prezentaci

Jaký je podle modelu s interakcí rozdíl zadní šířky samec − samice při délce krunýře 30 mm?

Tabulka čtyř koeficientů modelu s interakcí délky krunýře a pohlaví modrých krabů.

  • +1,97 mm
  • −0,12 mm
  • +5,71 mm
  • −1,76 mm

Správná odpověď: −1,76 mm

Koeficient pohlavisamec platí jen při délce 0 mm. Při 30 mm k němu přičteme 30 × interakční koeficient: 1,973 + 30 × (−0,1245) = −1,76 mm. Samec má tedy při 30 mm užší zadní část krunýře.

Jakou hodnotu má v modelu forma * pohlavi koeficient formaoranžová:pohlavisamec?

Průměrné délky krunýře samic a samců modré a oranžové formy; u modré formy jsou samci delší, u oranžové jsou obě pohlaví téměř stejně dlouhá.

  • −4,84 mm
  • −0,93 mm
  • +3,91 mm
  • +6,52 mm

Správná odpověď: −4,84 mm

Interakce je rozdíl dvou rozdílů: (33,69 − 34,62) − (32,01 − 28,10) = −0,93 − 3,91 = −4,84 mm. O tolik je rozdíl samec − samice u oranžové formy menší než u modré.

Co platí o přírůstku objemu na 1 cm průměru u stromu vysokého 35 m?

Tabulka čtyř koeficientů modelu objemu dřeva s interakcí průměru a výšky kmene a pozorované rozsahy obou prediktorů.

  • Je +0,107 m³ na 1 cm; vzorec platí pro každou výšku.
  • Je −0,065 m³ na 1 cm; to je koeficient prumer_cm.
  • Vyjde +0,107 m³ na 1 cm; 35 m je ale mimo data.
  • Nelze ho spočítat; u dvou číselných prediktorů sklon nemá smysl.

Správná odpověď: Vyjde +0,107 m³ na 1 cm; 35 m je ale mimo data.

Vzorec −0,065 + 35 × 0,00493 dá +0,107 m³ na 1 cm, ale změřené stromy byly vysoké jen 19–27 m. Mimo tento rozsah data tvar vztahu nepodporují, takže číslo je extrapolace, odhad mimo pozorované výšky. Uvnitř dat je sklon např. +0,038 (21 m) a +0,058 (25 m) m³ na 1 cm.

Výsledky učení

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

  • vypočítat R² z celkového součtu čtverců a součtu čtverců residuí a vysvětlit, co popisuje;
  • vysvětlit, proč nejvyšší R² automaticky neznamená nejlepší model;
  • rozlišit, jak R², adjustované R² a AIC zacházejí se shodou a složitostí, a používat AIC jen k relativnímu porovnání;
  • sestavit malou sadu biologicky zdůvodněných kandidátních modelů se stejnou odezvou a stejnými pozorováními;
  • obhájit volbu modelu pomocí biologické otázky, složitosti, diagnostiky a interpretovatelnosti a pojmenovat její omezení.

Jedna přímka pro všechny stromy?

Jeden řádek tabulky představuje jeden strom

data_stromy <-
  read.csv(file = "data/allometrie_stromu.csv")

data_stromy$druh_stromu <-
  factor(x = data_stromy$druh_stromu)

vec_promenne <-
  c("druh_stromu", "prumer_kmene_cm", "vyska_stromu_m")

head(
  x = data_stromy[, vec_promenne],
  n = 4
)
            druh_stromu prumer_kmene_cm vyska_stromu_m
1 Pseudotsuga menziesii           54.61          27.04
2 Pseudotsuga menziesii           34.80          27.42
3 Pseudotsuga menziesii           24.89          21.23
4 Pseudotsuga menziesii           28.70          24.96


63 stromů

  • Pinus monticola: 19 stromů
  • Pinus ponderosa: 22 stromů
  • Pseudotsuga menziesii: 22 stromů
  • chybějící hodnoty: 0
  • průměr kmene: 4,8–73,7 cm
  • výška: 3,6–45,0 m

Jiné stromy než minule: tři druhy jehličnanů a výška je tentokrát odezva.

Data: J. Marshall, University of Idaho · balíček lgrdata 0.1.1 (allometry) · CC0

Jak se výška mění s průměrem kmene?

Bodový graf výšky stromu proti průměru kmene pro jehličnany tří druhů; výška s průměrem roste.

Ve dvojici: Jaký vzorec v datech převládá? Co by mohlo vadit jedné přímce pro všechny stromy?

Jedna společná přímka pro všechny druhy


mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~
      prumer_kmene_cm,
    data = data_stromy
  )

coef(object = mod_spolecny)
    (Intercept) prumer_kmene_cm 
      7.5967386       0.5179077

Body stromů tří druhů a jedna společná fialová přímka modelu.

Sklon: o 0,52 m vyšší odhadovaná výška na každý centimetr průměru — stejně pro všechny druhy.

Co už umíme přečíst ze summary()?

summary(object = mod_spolecny)

Call:
lm(formula = vyska_stromu_m ~ prumer_kmene_cm, data = data_stromy)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.4226  -5.3191   0.7425   4.6637  10.5197 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      7.59674    1.47310   5.157 2.88e-06 ***
prumer_kmene_cm  0.51791    0.03647  14.201  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.548 on 61 degrees of freedom
Multiple R-squared:  0.7678,    Adjusted R-squared:  0.764 
F-statistic: 201.7 on 1 and 61 DF,  p-value: < 2.2e-16

Známe: volání modelu · residua · koeficienty s nejistotou

Zůstává v residuích rozdíl mezi druhy?

Residua společné přímky proti odhadnuté výšce, barevně podle druhu stromu.

Stejná residua s vodorovnými čarami průměrného residua každého druhu; Pinus monticola leží převážně nad nulou.


Leží residua některého druhu převážně nad nulou, nebo pod ní?

Průměrné residuum P. monticola: +2,7 m

Společná přímka přehlíží rozdíly mezi druhy.

Tři představy o růstu stromů

Jaké jiné přímky by mohly stromy popsat?

Stejný bodový graf stromů bez přímek, připravený pro náčrt alternativních modelů.


Ve dvojici: Načrtněte dvě jiné možnosti než jednu společnou přímku.

Co by každá z nich tvrdila o stromech?

Pomůcka: modely z minulé lekce — + a *.

Tři modely vyjadřují tři biologické představy

Tři schémata: jedna přímka, tři rovnoběžné přímky a tři přímky s různými sklony.

vyska ~ prumer

Všechny druhy sdílejí stejný intercept i sklon.

vyska ~ prumer + druh

Při stejném průměru se druhy mohou lišit výškou, sklon sdílejí.

vyska ~ prumer * druh

Vztah průměru a výšky se může mezi druhy lišit.

Tuto sadu kandidátů určíme předem — z biologické otázky, ne z výsledků.

Všechny tři modely odhadneme ze stejných 63 stromů

# Všechny tři kandidáty zapíšeme vedle sebe
mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm,
    data = data_stromy
  )

mod_aditivni <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm + druh_stromu,
    data = data_stromy
  )

mod_interakce <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm * druh_stromu,
    data = data_stromy
  )

Stejná odezva vyska_stromu_m · stejných 63 stromů v každém modelu

Který model prochází nejblíže pozorovaným stromům?

Tři panely se stejnými stromy: jedna společná přímka, tři rovnoběžné přímky a tři přímky s různými sklony.

A · Společná přímka

B · Posun mezi druhy

C · Různé sklony

Nejblíže procházejí přímky s různými sklony

Stejné tři panely se svislými oranžovými úsečkami residuí a součtem čtverců residuí v záhlaví; nejmenší je u modelu Různé sklony.

Je proto nejlepší?

Podle čeho poznáme lepší model?

Podle čeho byste mezi modely vybírali?

  • jak těsně model odpovídá pozorovaným datům;
  • kolik koeficientů potřebuje;
  • zda v residuích nezůstává systematický vzorec;
  • zda odpovídá biologické otázce;
  • zda jeho výsledek dokážeme smysluplně interpretovat.

Žádné jediné číslo neodpoví na všechny tyto otázky.

Jak změřit shodu modelu s daty?

Co už víme o součtu čtverců residuí?

Hledání přímky

Kandidátní přímky jsme porovnávali součtem čtverců residuí.

Zkratka RSS

Stejnou veličinu jsme označili RSS (residual sum of squares).

Skupinové průměry

Skupinové průměry RSS zmenšily oproti jednomu společnému průměru.


Co se stane s RSS, když přímka prochází blíže pozorovaným bodům?

Kratší residua → menší čtverce → menší RSS = těsnější shoda s daty.

Je přímka lepší než předpověď průměrem?

Šest skutečných stromů A až F a vodorovná přerušovaná čára jejich průměrné výšky.


Šest skutečných stromů z tabulky.

Předpověď bez prediktoru: každému stromu průměrnou výšku.


Prošla by přímka podle průměru kmene body blíže než průměr?

Jak velké jsou rozdíly výšek kolem průměru?

Stejných šest stromů se svislými úsečkami od průměrné výšky k pozorované výšce.

Strom Odchylka Čtverec
A 1,7 2,9
B 6,5 42,0
C −6,0 35,5
D −2,0 4,0
E −19,5 380,2
F 19,3 372,2

Součet čtverců: 836,8 m²

Tomuto součtu říkáme celkový součet čtverců (total sum of squares).

Kolik rozdílů zůstane po odhadnutí přímky?

Stejných šest stromů s fialovou odhadnutou přímkou a oranžovými úsečkami residuí.

Strom Residuum Čtverec
A −6,6 43,9
B 3,6 12,8
C 1,7 3,0
D 6,4 41,4
E −6,6 44,1
F 1,5 2,3

RSS: 147,6 m² · celkový: 836,8 m²

Který strom leží nejvýše nad přímkou? A který nejvíce přispívá k RSS? Strom D; k RSS nejvíce E (rozhoduje velikost, ne znaménko).

Jakou část rozdílů přímka zachytila?

\[ \text{zachycený podíl} = 1 - \frac{\text{součet čtverců residuí}}{\text{celkový součet čtverců}} \]

Jaká část rozdílů ve výšce zůstala a jakou část přímka zachytila? Dosaďte RSS 147,6 m² a celkový součet čtverců 836,8 m².

\[ \text{zachycený podíl} = 1 - \frac{\text{147,6}}{\text{836,8}} \]

\[ \text{zachycený podíl} = \text{0,824} \]

Zachycený podíl se jmenuje R²

Koeficient determinace R² (coefficient of determination)

\[ \Large R^2 = 1 - \frac{SS_{\mathrm{res}}}{SS_{\mathrm{total}}} \]

  • \(SS_{\mathrm{res}}\) — součet čtverců residuí (RSS) od odhadnuté přímky;
  • \(SS_{\mathrm{total}}\) — celkový součet čtverců odchylek od průměru;
  • u šesti stromů: \(R^2 = \text{0,824}\).

Jak se mění RSS a R², když přímku natáčíme?

Animace: šest stromů a fialová přímka, která se natáčí kolem středu dat od vodorovného průměru k odhadnuté přímce a zpět. Oranžové úsečky residuí se zkracují, sloupec RSS se zmenšuje vůči celkovému součtu čtverců a R² roste od 0 k hodnotě odhadnuté přímky.

Sledujte oranžové úsečky a sloupec RSS: kdy je R² nulové a kdy největší?

Stejných šest stromů, čtyři natočení přímky

Čtyři panely se stejnými šesti stromy: přímka natočená na 0, třetinu, dvě třetiny a celý odhadnutý sklon. Residua se zkracují, RSS klesá a R² roste od 0 k hodnotě odhadnuté přímky.

Vodorovná přímka = předpověď průměrem: R² = 0. Odhadnutá přímka: nejmenší RSS, a proto největší R².

Dává summary() stejné R² pro všech 63 stromů?

\[ R^2 = 1 - \frac{\text{1877,4}}{\text{8084,6}} = \text{0,768} \]

summary(object = mod_spolecny)$r.squared
[1] 0.7677804

Multiple R-squared ve výpisu je tentýž výpočet — nic nového.

Společná přímka má R² = 0,768. Co to znamená?

A

Model správně předpověděl 76,8 % stromů.

B

Model zachytil 76,8 % variability výšky kolem průměru v těchto datech.

C

S pravděpodobností 0,768 je model pravdivý.

Přesvědčte souseda, proč jste vybrali svoji možnost.

Obhajitelná je B: R² popisuje podíl zachycené variability v datech, ze kterých byl model odhadnut.

Co R² popisuje a na co nestačí?

R² popisuje

  • podíl variability odezvy zachycený modelem;
  • ve stejných datech, na kterých byl model odhadnut.

R² neříká

  • zda model splňuje předpoklady;
  • zda je vztah příčinný;
  • jak dobře model poslouží pro nové stromy;
  • zda je to nejlepší z našich kandidátů?

Pomáhá každý koeficient navíc?

Co udělá R² po přidání druhu a interakce?

prumer → prumer + druh → prumer * druh

Co se stane s R², když přidáme druh a interakci?

A · R² poroste

B · R² může klesnout

C · R² zůstane stejné

R² po přidání členů neklesne

R² tří kandidátních modelů roste od společné přímky přes posun mezi druhy k různým sklonům.


Model k RSS
Společná přímka 2 1877
Posun mezi druhy 4 1675
Různé sklony 6 1644


Složitější model umí zopakovat jednodušší, a tak R² nejhůř zůstane stejné; klesnout nemůže.

Vyšší R² není automatické vítězství. I téměř zbytečný prediktor R² nesníží.

Které koeficienty jednotlivé modely odhadují?

Který seznam patří společné přímce, který posunu mezi druhy a který různým sklonům?

A· Posun mezi druhy · k = 4

(Intercept)
prumer_kmene_cm
druh_stromuPinus ponderosa
druh_stromuPseudotsuga menziesii

B· Společná přímka · k = 2

(Intercept)
prumer_kmene_cm

C· Různé sklony · k = 6

(Intercept)
prumer_kmene_cm
druh_stromuPinus ponderosa
druh_stromuPseudotsuga menziesii
prumer_kmene_cm:druh_stromuPinus ponderosa
prumer_kmene_cm:druh_stromuPseudotsuga menziesii

Jak R² opravit za počet koeficientů?

\[ \text{adjustované R}^2 = 1 - \text{nezachycený podíl} \times \text{úprava za počet stromů a koeficientů} \]

Posun mezi druhy: \(R^2 = \text{0,793}\), \(n = 63\) stromů, \(k = 4\) koeficienty

\[ R^2_{\mathrm{adj}} = 1 - \left(1 - \text{0,793}\right) \times \frac{63 - 1}{63 - 4} = \text{0,782} \]

Adjustované R² (adjusted R-squared):

\[ \Large R^2_{\mathrm{adj}} = 1 - \left(1 - R^2\right) \times \frac{n - 1}{n - k} \]

Co udělá adjustované R² po přidání interakce?

Posun mezi druhy → Různé sklony: R² 0,793 → 0,797 · k 4 → 6 · n = 63

A · adjustované R² vzroste

B · adjustované R² klesne

Adjustované R² po přidání interakce klesne

R² tří modelů stále roste, zatímco adjustované R² je nejvyšší u posunu mezi druhy a u různých sklonů klesne.


Dva sklony navíc zlepšily R² jen o 0,004.


Zlomek \(\frac{n-1}{n-k}\) vzrostl víc, než se zmenšil nezachycený podíl → úprava zvítězila.

Co už o třech modelech víme?

Jak poznáme, že je model dobrý?

Už umíme říct

  • R² měří shodu s daty a vždy upřednostní nejpružnější model.
  • Adjustované R² upřednostní posun mezi druhy.

Zatím nevíme

  • Kdy se složitost obecně vyplatí?
  • Vypadají residua vybraného modelu v pořádku?

Kdy se složitost vyplatí?

Vyplatí se model šitý na míru?

Liška krejčí šije jehličnanu oblek s desítkami švů a špendlíků, které kopírují každý hrbol kmene; vedle visí jednoduchý dobře padnoucí kabát a za nimi čeká další strom.

Ilustrace vytvořená pomocí AI.

Ve dvojici: Posun mezi druhy a různé sklony popisují stromy téměř stejně dobře (R² 0,793 a 0,797). První potřebuje 4 koeficienty, druhý 6.

Který byste upřednostnili? Co byste chtěli vědět, než se rozhodnete?

Princip úspornosti (parsimony)

Složitější model upřednostníme jen tehdy, když dostatečně zlepší shodu s daty.

AIC váží shodu i složitost

Akaikeho informační kritérium (Akaike information criterion, AIC)

  • lepší shoda AIC snižuje;
  • každá odhadovaná veličina navíc AIC zvyšuje;
  • preferujeme nižší AIC.
AIC(object = mod_spolecny)
[1] 398.64

AIC není podíl a nemá stejnou škálu jako R².

Model má AIC 398,64. Je to dobrý model?

A · Ano, AIC je nízké.

B · Ne, AIC je vysoké.

C · Z jedné hodnoty to nepoznáme.

AIC je relativní. Porovnáváme jen modely se stejnou odezvou, stejnými pozorováními a stejným předpokladem o rozdělení residuí.

O kolik je AIC každého kandidáta vyšší než nejnižší?

\[ \Delta\text{AIC} = \text{AIC kandidáta} - \text{nejnižší AIC v sadě} \]

Společná přímka: \(\text{398,64} - \text{395,44} = \text{3,20}\)

Model AIC ΔAIC
Společná přímka 398,64 3,20
Posun mezi druhy 395,44 0,00
Různé sklony 398,27 2,83

Posun mezi druhy má ΔAIC = 0 · univerzální práh nepoužíváme

Ukazují R², adjustované R² a AIC na stejný model?

Model k R² ↑ adj. R² ↑ AIC ↓
Společná přímka 2 0,768 0,764 398,64
Posun mezi druhy 4 0,793 0,782 395,44
Různé sklony 6 0,797 0,779 398,27


Ve dvojici: Proč R² upřednostní jiný model než adjustované R² a AIC?

R² hodnotí jen shodu. Adjustované R² a AIC penalizují složitost — dva sklony navíc se nevyplatily. Nejsou to tři hlasy o pravdě.

Nejlépe podporovaný kandidát nemusí být dobrý model

Porota lesních zvířat s prázdnými hodnoticími kartičkami posuzuje na pódiu tři kandidáty, každého s jiným uspořádáním rovných tyček; za pootevřenými dveřmi čeká čtvrtý kandidát s prohnutou tyčkou.

Ilustrace vytvořená pomocí AI.

Porovnání řadí jen modely, které jsme do sady zařadili.

Neříká, že je vítěz obecně dobrý. Model, který jsme vůbec neposuzovali, může stromy popsat lépe.

Proto vybraný model ještě zkontrolujeme.

Obstojí vybraný model?

Zůstává v residuích modelu s posunem mezi druhy nějaký vzorec?

Vlevo residua modelu s posunem mezi druhy proti odhadnuté výšce, vpravo Q–Q graf residuí.

Stejné grafy; v residuích je navíc vyhlazená oranžová křivka, která tvoří oblouk: nízké i nejvyšší odhadnuté výšky mají residua převážně pod nulou, středně vysoké nad nulou.

Ve dvojici: Leží residua v celém rozsahu odhadnutých výšek náhodně kolem nuly?

Oranžová křivka jen vyhlazuje průběh bodů, aby byl vzor lépe vidět: residua tvoří oblouk · Q–Q graf není přímý.

Které stromy ovlivňují model nejvíc?

Cookovy vzdálenosti všech stromů pro model Posun mezi druhy s přerušovanou orientační čarou 4/n; několik stromů leží nad ní.

Cookova vzdálenost: o kolik by se odhadnuté hodnoty změnily, kdyby strom chyběl. Nad orientační čarou \(4/n\): 4 stromy.

Důvod stromy prověřit — ne je automaticky odstranit.

Obhajte model dvěma větami

40 s sami · ve dvojici porovnejte

„Ze tří kandidátů vybíráme model ____, protože ____.“

„Omezením tohoto závěru je ____.“

Opřete se aspoň o dva různé druhy důkazů: shodu, složitost, diagnostiku, biologický smysl.

Který model lze pro tyto stromy nejlépe obhájit?

Posun mezi druhy

  • zachytí rozdíly mezi druhy, které společná přímka přehlíží;
  • nejvyšší adjustované R² (0,782) a nejnižší AIC.

Omezení

  • vybírali jsme jen ze tří kandidátů;
  • residua tvoří oblouk: přímka nevystihuje tvar vztahu;
  • 63 stromů, Q–Q graf není přímý;
  • pozorovací data: nejde o příčinný efekt druhu.

Dobrý model nepoznáme podle jediného čísla. Obhajoba propojuje otázku, shodu, složitost, diagnostiku a interpretovatelnost.

Jak porovnávat modely bez automatického hledání?

  1. Položte biologickou otázku.
  2. Prohlédněte a zobrazte data.
  3. Určete malou sadu zdůvodněných kandidátů.
  4. Použijte stejnou odezvu a stejná pozorování.
  5. Porovnejte přímky, residua, R², adjustované R² a AIC.
  6. Zkontrolujte předpoklady a vlivná pozorování.
  7. Obhajte závěr a pojmenujte omezení.


Automatické přidávání a odebírání členů vytvoří mnoho modelů, ale neurčí, které dávají biologický smysl.

Co si odnést

  • R² = podíl zachycené variability; počítá se z celkového součtu čtverců a RSS.
  • Po přidání členu R² neklesne, a proto samo model nevybírá.
  • Adjustované R² a AIC penalizují složitost, každé jinak; AIC je relativní.
  • Kandidáty určujeme předem: stejná odezva, stejná pozorování.
  • Obhajoba spojuje několik druhů důkazů včetně diagnostiky a omezení.

Data: J. Marshall, University of Idaho · lgrdata::allometry · CC0

Co když žádná přímka nestačí?

Řada jehličnanů od mladých po staré: jejich výška nejprve roste a pak se ustálí, zatímco obří rovná fialová tyč stoupá dál až do mraků; liška a veverka se na ni udiveně dívají.

Ilustrace vytvořená pomocí AI.

Všechny tři kandidáty jsme kreslili jako přímky — a residua přesto tvoří oblouk.

Jak popsat vztah, který není přímkový?