Lineární modely (stručné shrnutí) (Obecné) lineární modely • Jedna odpověď, jeden a více prediktorů - Prediktory kvantitativní i kategoriální - ANOVAs prediktorem o n úrovních je analogická mnohonásobné regresi o n-1 prediktorech • Rovnice: y = a + bx1 + cx2 + ...+£ - y: odpověď - a: intercept - b,c: regresní koeficienty prediktorů - z: reziduály (residua): předpokládá se, že všechny pocházejí ze stejného normálního rozložení N(0, a) Modely s více prediktory • Dvoucestná (Vícecestná) ANOVA - odpověď ~ faktor. 1 + faktor.2 + ... • Mnohonásobná regrese - odpověď - lin.prediktor.1 + lin.prediktor.2 + ... • Lineárni model (Analýza kovariance): kategoriální i lineárni prediktory • Aditivní efekty prediktorů vs. interakce - odpověď na faktor. 1 (prediktor.1) závisí na hodnotě faktoru.2 (prediktoru.2) - aditivitu lze statisticky testovat a případně zamítnout ve prospěch interakce Interakce mezi prediktory •Průkazná interakce znamená vzájemné ovlivňování vlivu prediktorů na odpověď -y = a + bx1 + cx2 + dx1x2 + £ •Test interakce HO: d = 0 -d > 0: pozitivní int., vyšší hodnoty než aditivita -d < 0 : negativní int., nižší hodnoty než aditivita •dfint = dfx1*dfx2 •Jak napsat -x (Alt + 0215) formálne v odborných textech -: v R pouze interakční člen -*vR znamená additivitu dohromady s interakcí -rovnici nahoře v R napíšeme takto y ~ x1*x2 •Interaction plot Znázornění interakce -Pokud int. není přítomna, budou čáry paralelní Interakce neznamená korelaci prediktorů! —i LU LU m E o to Cl Interaction plot o o CN O O O O Li") O - Control x.factor data biomass Reakce růstu rostlin na zalévání a hnojení Experimentální design: 24 květináčů s rostlinami rozdělenými náhodně do 4 skupin s faktoriálním uspořádáním hnojení a velikosti zálivky Otázka: Jaký je vliv hnojení a zalévání na produkci nadzemní biomasy rostlin? aov.l<-aov(biomass-watering*fertil, data=plants) summary(aov.1) watering fertil 65 Control Control 58 Control Control 74 Control Control 65 Control Control 81 Control Control 78 Control Control 92 Watered Control 86 Watered Control 94 Watered Control 100 Watered Control 89 Watered Control 90 Watered Control 110 Control Fertilized 118 Control Fertilized 128 Control Fertilized 121 Control Fertilized 119 Control Fertilized 116 Control Fertilized 185 Watered Fertilized 196 Watered Fertilized 201 Watered Fertilized 195 Watered Fertilized 193 Watered Fertilized 191 Watered Fertilized Df Sum Sq Mean Sq F value Pr(>F) watering 1 13968 13968 336.38 5.604e- ]_ 4 * * * fertil 1 33825 33825 814.57 < 2.2e- ]_ g -k-k-k watering:fertil 1 4240 4240 102.11 2.654e- Q g -k -k -k Residuals 20 831 42 Interaction plot Cl Control Fertilized x.f a cto r data biomass Reakce růstu rostlin na zalévání a hnojení Experimentální design: 24 květináčů s rostlinami rozdělenými náhodně do 4 skupin s faktoriálním uspořídáním hnojení a velikosti zálivky Otázka: Jaký je vliv hnojení a zalévání na produkci nadzemní biomasy rostlin? aov.l<-aov(biomass-watering*fertil, data=plants) summary(aov.1) Df Sum Sq Mean Sq F value Pr(>F) 336.38 5.604e-14 *** 814.57 < 2.2e-16 *** 102.11 2.654e-09 *** Interaction plot Závěr: Hnojení i zalévání mají průkazný pozitivní vliv růst rostlin. Mezi prediktory je navíc pozitivní interakce. watering fertil 65 Control Control 58 Control Control 74 Control Control 65 Control Control 81 Control Control 78 Control Control 92 Watered Control 86 Watered Control 94 Watered Control 100 Watered Control 89 Watered Control 90 Watered Control 110 Control Fertilized 118 Control Fertilized 128 Control Fertilized 121 Control Fertilized 119 Control Fertilized 116 Control Fertilized 185 Watered Fertilized 196 Watered Fertilized 201 Watered Fertilized 195 Watered Fertilized 193 Watered Fertilized 191 Watered Fertilized watering 1 13968 13968 fertil 1 33825 33825 watering:fertil 1 4240 4240 Residuals 20 831 42 E ^ Q_ Control Fertilized x.f actor Testování vlivu jednotlivých prediktorů / výběr prediktorů do modelu • Cíle - ukázat, které prediktory mají průkazný vliv • lze testovat jednoduché (marginální) efekty, tj. párové korelace mezi jednotlivými prediktory a odpovědí - sestavit minimální adekvátní model, který bude zahrnovat pouze průkazné prediktory • kondicionální (parciální) efekty prediktorů • Statistické testy tohle moc neumí - umí otestovat model, porovnat kvalitu více modelů • test signifikance • AIC Akaike information Criterion • Kvantifikuje množství informaci v odpovědi, kterou vysvětluje model - umožňuje porovnání dvou modelů, které se liší počtem stupňu volnosti (modelu) - nižší AIC značí lepší fit než vyšší AIC (absolutní hodnoty nejsou důležité) • AIC = 2k - 2log(L), (log = přirozený log), k je počet parametrů (tedy df model v lin. modelu) - v lin. mod. AIC = 2k - 2 log (n/RSS) + C, kde RSS je reziduálni suma čtverců, C je konstanta (lze ignorovat) • Různé názory na možnost kombinovat s F-testem průkaznosti Výběr prediktorů do modelu • Postupný - Forward selection: k nulovému modelu přidávám postupně prediktory • vhodnější pro observační data - Backward selection: ze saturovaného modelu odebírám nevýznamné prediktory • vhodnější pro experimentální data - Oba směry: zvažují v každém kroku přidání I odebrání prediktorů Problém s mnohonásobným paralelním testováním • Provedu-li více testů s pravděpodobností chyby I. druhu 0.05, pravděpodobnost, že udělám chybu aspoň jednou velmi vzrůstá. - pro dva testy p = 0.05 + 0.05 - 0.05 x 0.05 = 0.0975 • Řešení: různé korekce (Holm, Bonferroni, false detection rate -FDR) upravují p-hodnoty nahom, čímž kontrolují/redukují riziko chyby • Alternativa: "protected multiple testing" - spočtu test saturovaného modelu se všemi prediktory. Je-li průkazný, další korekci kvůli multiple testing už neřeším - Je-li ovšem neprůkazný, tak s testováním končím se závěrem, že regrese není průkazná. • Skutečný problém hlavně v analýzách "velkých" dat z databází. Lineární model pro experimentální data Sample seedlings 1 2 3 4 5 6 7 8 9 10 11 treatment 7 control 5 control 12 control 8 control 13 control 12 control 7 control 5 control 7 control 6 control 19 grazing productivity temperature 714 7.2 518 4.5 379 7.4 686 5.4 703 5 775 6 651 6.2 630 7.6 470 8.4 557 4.7 394 6.8 to be continued Q: Jaký je vliv obhospodařování louky (kosení, pastva) na počet semenáčků, které se objeví na jaře? Ovlivňují semenáčky ještě další faktory prostředí? Design: 30 experimentálních ploch (10 na každý typ obhospodařování) náhodně rozmístěné v krajině (různé lokality). Control = opuštěná nekosená louka. Zaznamenána byla i produktivita a průměrná teplota lokality. lm.full<-lm(seedlings-treatment*productivity*temperature, data=seedl] anova(lm.full) # Analysis of Variance Table # Response: seedlings # # # # # # # # # # # # Df treatment 2 productivity 1 temperature 1 treatment:productivity 2 treatment:temperature 2 productivity:temperature 1 treatment:productivity:temperature 2 Residuals 18 Signif. codes 0 0 . 001 Sum Sq Mean Sq F value 618.07 309.033 31.5710 290.61 290.605 29.6883 1.88 1.882 0.1922 5.02 2.509 0.2563 6.87 3.434 0.3508 0.25 0.250 0.0256 6.48 3.242 0.3312 176.19 9.789 0.01 0.05 0.1 * Pr (>F) , 301e-06 , 553e-05 ~k ~k ~k ~k ~k ~k 0 0 0 0 0 6663 7767 7088 8747 7223 Tuto tabulku lze použít jako výstup testu jednotlivých prediktorů a jejich interakcí Lineární model pro experimentální data Další krok: odstraňte neprůkazné prediktory z modelu (pomocí backward selekce; ponechejte I neprůkazné main efekty pokud je průkazná interakce) lm.final<-lm(seedlings-treatment+productivity, data=seedl) summary(lm.final) Call: lm(formula = seedlings ~ treatment + productivity, data = seedl) Residuals Min -4.698 -1 1Q Median !40 -0.315 3Q 1. 975 Max 4 .741 Coefficients: (Intercept) treatmentgrazing treatmentmowing productivity Signif. codes: 0 Estimate Std. Error t value Pr(>|t|) 22.315470 2.437909 9.154 1.29e-09 *** 9.581024 1.251180 7.658 3.98e-08 *** 2.173362 1.268726 1.713 0.0986 . -0.024764 0.003996 -6.198 1.48e-06 *** y -k-k-k r 0.001 y -k-k r 0.01 y -k r 0 . 05 0.1 Residual standard error: 2.75 on 26 degrees of freedom Multiple R-squared: 0.8221, Adjusted R-squared: 0.8015 F-statistic: 40.04 on 3 and 26 DF, p-value: 6.857e-10 Koeficienty prediktorů (nebo kontrasty pro faktory) Celkový test modelu Lineární model pro experimentální data Další krok: odstraňte neprůkazné prediktory z modelu (pomocí backward selekce; ponechejte I neprůkazné main efekty pokud je průkazná interakce) lm.final<-lm(seedlings-treatment+productivity, data=seedl) summary(lm.final) Call: lm(formula seedlings ~ treatment + productivity, data = seedl Residuals Min -4.698 -1 1Q Median :40 -0.315 3Q 1. 975 Max 4 .741 Coefficients: (Intercept) treatmentgrazing treatmentmowing productivity Estimate Std. Error t value Pr(>|t|) 22.315470 2.437909 9.154 1.29e-09 *** 9.581024 1.251180 7.658 3.98e-08 *** 2.173362 1.268726 1.713 0.0986 . -0.024764 0.003996 -6.198 1.48e-06 *** Koeficienty prediktorů (nebo kontrasty pro faktory) Signif. codes: 0 0.001 0.01 0.05 0.1 Residual standard error: 2.75 on 26 degrees of freedom Multiple R-sguared: 0.8221, Adjusted R-sguared: 0.8015 F-statistic: 40.04 on 3 and 26 DF, p-value: 6.857e-10 Celkový test modelu Závěr: Průkazný vliv na počet semenáčků má způsob obhospodařování a produktivita mají. Jejich efekty jsou aditivní. Produktivita snižuje počet semenáčků. Obhospodařivání jejich počet zvyšuje. Průkazné zvýšení je ale způsobené pouze pastvou. Lineární model pro observační data • Typicky mnoho potenciálních prediktorů - mnohdy nemožnost fitovat saturovaný model (potenciálních prediktorů více než pozorování) • Forward selection s korekcí - akceptovatelná varianta • Korelované potenciální prediktory: těžko řešitelný problém - např. v ČR: nadmořská výška, průměrná teplota a úhrn srážek - Možnost prezentovat jednoduché korelace spolu s finální lineárním modelem Výběr prediktorů do lineárního modelu • Spousta možností jak to udělat v zásadě dobře - postupný výběr s korekcí (ať už jakoukoliv) - využití "test protection" • Jednoznačně špatně - Statistical fishing - Forward selekce bez korekce pro observační data • speciálně při použití mnoha potenciálních prediktorů - Úpravy testů "aby to vyšlo průkazně"