Stochastické modely časových řad (2. cvičení) 1 M5201 - 2. CVIČENÍ: Regresní analýza v R 1 Úvod do regresní analýzy Regresní analýza je jednou z nejčastěji používaných statistických technik. Zabýváme se při ní vztahy mezi několika veličinami. Snažíme se zjistit, jak hodnoty jedné veličiny (vysvětlované proměnné, odezvy) závisí na hodnotách ostatních veličin (vysvětlujících veličin, regresorů, kovariát). Vycházíme z předpokladu, že vysvětlovaná veličina je funkcí vysvětlujících proměnných a náhodné složky, což je náhodná veličina s nulovou střední hodnotou, která reprezentuje náhodné odchylky, chyby v měření apod. Y = f(x1,x2,.. .,xm,s), kde Y je vysvětlovaná proměnná xi,..., xm jsou vysvětlující proměnné e je náhodná složka Funkce / patří do předem zvolené třídy funkcí a závisí na jednom či více parametrech. Ty jsou neznámé a odhadují se na základě zjištěných dat. Obvykle máme k dispozici n pozorování vysvětlované proměnné Y\, Y2,..., Yn a jim odpovídající hodnoty vysvětlujících proměnných xn,x12,xln, x2i,x22, x2n, %i, • • •, xmn a pomocí modelu můžeme popsat jednotlivá pozorování Y{ — f(x\i, X2{, • • •, xmi, £j), i — 1,..., n Nejjednodušším případem regresního modelu je lineární regresní model. Je charakteristický tím, že vysvětlovanou proměnnou se v něm snažíme popsat pomocí lineární kombinace známých funkcí vysvětlovaných proměnných s neznámými regresními koeficienty. Pomocí rovnice můžeme lineární regresní model zapsat takto Yi = fa + f3ifi(xu,x2i,..., xmi) + f32f2(xu,x2i,..., xmi) H----Pkík(xu-,x2i,..., xmi) + a, kde Y je vysvětlovaná proměnná xi,..., xm jsou vysvětlující proměnné /i) f2, ■ ■ ■, f k jsou známé (předem zvolené) funkce vysvětlujících proměnných 00, f3i,..., (3k jsou neznámé regresní koeficienty £ je náhodná veličina s nulovou střední hodnotou a konstatním rozptylem V lineárním regresním modelu tedy vysvětlovanou proměnnou modelujeme jako součet systematické složky (což je funkce, která je lineární v neznámých regresních koeficientech) a náhodné složky. V lineárním regresním modelu se dále předpokládá, že pro různé realizace Yi, i = 1,... ,n vysvětlované proměnné jsou odpovídající náhodné veličiny eí navzájem nezávislé a mají stejné rozdělení s nulovou střední hodnotou a konstatním rozptylem, což značíme e j ~ IID(0, a2). Pro odhad neznámých parametrů regresních modelů existuje více metod. Pro lineární regresní modely se obvykle používá metoda nejmenších čtverců. Odhady parametrů se obvykle 2 Stochastické modely časových řad (2. cvičení) označují stříškou, např. /3o,/?i,.... Hodnoty vysvětlované proměnné Yi, které získáme po dosazení odhadnutých parametrů do regresní rovnice Yi = A) + /3ifi{xii,x2i,xmi) + f32f2(xu,x2i,..., xmi) H----l3kfk{xu,x2i,..., xmi) se nazývají vyrovnané hodnoty. Rozdíly mezi skutečnými a vyrovnanými hodnotami vysvětlované proměnné e% = Yi — Yi se nazývají rezidua. Veličiny v modelu mohou být jak spojité, tak kategorické (sem patří nominální, ordinální a diskrétní veličiny). Kategorické proměnné se také nazývají faktory a jejich hodnoty se nazývají úrovně. Značení: Spojité veličiny budeme značit písmeny z konce abecedy: x, Y,... Faktory budeme značit písmeny ze začátku abecedy: A,B,... 2 Použití funkce lm K odhadu parametrů lineárního regresního modelu je v systému R k dispozici funkce lm. Tato funkce má následující argumenty lm(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALŠE, y = FALŠE, qr = TRUE, singulár.ok = TRUE, contrasts = NULL, offset, ...) Nejdůležitější je argument formula, kterým se zadává, jakou podobu má mít rovnice popisující model. Modelový vztah mezi proměnnými se obecně zadává ve tvaru vysvětlovaná proměnná ~ systematická složka Zatímco na levé straně můžeme normálně používat obvyklé matematické funkce jako logaritmus nebo matematické operátory +, — apod., na pravé straně mají matematické operátory poněkud odlišný význam. Konkrétně: + přidat vysvětlující proměnnou odstranit vysvětlující proměnnou : interakce mezi veličinami * všechny komponenty 1 absolutní člen -1 odebrat absolutní člen Pokud potřebujeme při specifikaci modelu použít matematické operátory v obvyklém smyslu, musíme výraz s těmito operátory napsat jako argument funkce I(). Rozdíl můžeme ukázat na příkladu. y ~ xl + x2 definuje model Yi = /3q + PiXu + /32x2i + eí, zatímco y ~ l(xl + x2) definuje model Yi = /3q + /3i(xu + x2i) + £j. Konkrétní příklady zadávání modelů v systému R jsou ukázány v následující tabulce. Stochastické modely časových řad (2. cvičení) 3 Formální zápis Syntaxe v R Popis Yt = fi + St y ~ 1 Model pouze s absolutním členem Yi = Po + PiXi + el y ~ x Lineární model se spojitou vysvětlující proměnnou x (regresní přímka) y ~ x - 1 Lineární model se spojitou vysvětlující x bez absolutního členu log(Yi) = Po + PiXi + el log (y) " x Lineární model se spojitou vysvětlující x a s logaritmicky transformovanou vysvětlovanou proměnnou Yi = Po + PiXi + faxf + £i y ~ x + I(x' -2), Kvadratický model se spojitou y ~ poly(x, 2) vysvětlující proměnnou x (regresní parabola) Yi = Po + Pixu + p2x21 + £t y ~ xl+x2 Model se dvěma spojitými vysvětlujícími proměnnými x\ a x2, v každé z nich je lineární (vícenásobná regrese) Yi = Po + Pixu + p2x21 + y ~ xl*x2, Model se dvěma spojitými Pi,xllx2l + et y ~ xl+x2+I(xl*x2) vysvětlujícími proměnnými x\ a x2 s křížovým členem Yij — n + at + Bij y ~ A Model s jedním faktorem (jednofak-torová analýza rozptylu, ANOVA) Yijk — fJ, + ai + Pj + eljk y ~ A + B Model se dvěma faktory bez interakce (dvoufaktorová analýza rozptylu) Yljk = iJ,+al+Pj + (aP)lj+eljk y ~ y ~ A + B + A * B A:B, Model se dvěma faktory s interakcí (dvoufaktorová analýza rozptylu) Yijki — n + ai+ Pj + 7fc + Eijki y ~ A + B + C Model se třemi faktory bez interakcí (třífaktorová analýza rozptylu) Yijki = p, + cti + Pj + 7fc + y ~ A+B+C+A: :B+A: :C+B:C+A: : B: C Model se třemi faktory s in- (a/% + («7)* + (Pl)jk + y ~ A*B*C terakcemi (třífaktorová analýza {aP^)ljk + Eijki rozptylu) Yijki = M + cti + Pj + lk + y ~ A+B+C+A: :B+A: :C+B:C Model se třemi faktory obsahu- (aP)lj + (a-f)lk + (P'f)jk+eljki y ~ (A+B+C)' jící kromě hlavních efektů maximálně dvojné interakce (třífaktorová analýza rozptylu) Yij — H + Oj + fiXij + Bij y ~ x+A Model se spojitou proměnnou x a faktorem A bez interakce interakce (ANCOVA) — H~ j ~h j3xij -\- ôj Xij ~h S í j y ~ x*A Model se spojitou proměnnou x a faktorem A s interakcí (ANCOVA) Tabulka 1: Přehled zápisu základních modelů 4 Stochastické modely časových řad (2. cvičení) 3 Příklady Příklad 1: Regresní přímka Máme k dispozici údaje o tom, jaký počet australských domácností měl v letech 1998-2008 připojení na internet: Rok 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 Počet domácností 1098 1538 2340 3114 3445 4039 4393 4730 5138 5492 5878 Data nejprve načteme, a to tak, že údaje o počtech domácností zkopírujeme do schránky a použijeme funkci scan. > domácnosti <- scan("clipboard", sep = " ") > time <- 1998:2008 Načtená data vykreslíme do grafu. > TxtX <- "Rok" > TxtY <- "Počet domácnosti" > plot (time, domácnosti, type = "p", xlab = TxtX, ylab = TxtY) Q- o n-1-1-1-1-r 1998 2000 2002 2004 2006 2008 Rok Obrázek 1: Počet domácností s připojením na internet v letech 1998-2008 Budeme uvažovat regresní model, v němž bude počet domácností majících připojení na internet záviset lineárně na čase. Model můžeme popsat rovnicí Yi = fa + PiXi + Si 1, - - -, 11 kde Y{ je počet domácností s připojením na internet v i-tém roce Xi je rok odpovídající i-tému pozorování A)> Pí jsou neznámé regresní koeficienty, které budeme odhadovat Ei je náhodná odchylka příslušející i-tému pozorování Parametry tohoto modelu odhadneme v systému R pomocí funkce lm a k zobrazení výsledků použijeme funkci summary. Stochastické modely časových řad (2. cvičení) 5 > modeli <- lm(domacnosti ~ time) > summary (modeli) Call: lm(formula = domácnosti ~ time) Residuals: Min 1Q Median 3Q Max -306.45 -200.05 20.18 173.09 318.82 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -948407.45 45085.31 -21.04 5.81e-09 *** time 475.36 22.51 21.12 5.61e-09 *** Signif. codes: 0 "***' 0.001 "**' 0.01 "*' 0.05 ".' 0.1 x ' 1 Residual standard error: 236.1 on 9 degrees of freedom Multiple R-squared: 0.9802, Adjusted R-squared: 0.978 F-statistic: 446 on 1 and 9 DF, p-value: 5.615e-09 Ve výpisů výsledků odhadu modelu vidíme nejprve zadaný tvar modelu. Dále jsou zde uvedeny informace o reziduích (minimální a maximální reziduum, horní a dolní kvartil a medián). V další tabulce jsou shrnuty odhadnuté parametry a v dalších sloupcích jsou směrodatné odchylky odhadnutých parametrů a testové statistiky pro testování některých hypotéz týkajících se parametrů modelu. V závěru jsou uvedeny souhrnné statistiky týkající se celého modelu. Na závěr zobrazíme data proložená odhadnutou regresní přímkou. > plot(time, domácnosti, type = "p", xlab = TxtX, ylab = TxtY) > abline(lm(domacnosti ~ time), col = "blue", lty = 1) A) 948407.45 (Intercept) pi = 475.36 (time) o o o CD O o o o o 1998 2000 2002 2004 2006 2008 Rok Obrázek 2: Počet domácností s připojením na internet v letech 1998-2008 6 Stochastické modely časových řad (2. cvičení) příklad 2: polynomická regrese Máme k dispozici průměrné roční průtoky vody v řece Nigeru v Coulicouro (Mali). Údaje se vztahují k období 1907 až 1957 a jsou uvedena v jednotkách cfs.l0~3. Data jsou uložena v souboru DataNiger.dat. První řádek obsahuje popis časové řady, pak následuje volný řádek, teprve potom dva sloupce dat: rok a průměrný roční průtok. Nejprve načteme nadpis. > íileDat <- paste (data. library, "DataNiger.dat", sep = "") > (TXT <- paste(scan(fileDat, what = "", nlines = 1), collapse = " ")) [1] "Průměrné rocni průtoky vody v~rece Nigeru v~Coulicoure (Mali) v~letech 1907 az 1957" Pak načteme vlastní data do datového rámce za pomoci funkce read.table s argumentem skip=2, který způsobí vynechání prvních dvou řádků. Ve vzniklém datovém rámci pojmenujeme proměnné, vypíšeme jeho strukturu a prvních šest řádků. Nakonec data vykreslíme. > dataNiger <- read.table(ÍileDat, header = F, skip = 2) > names (dataNiger) <- c ("Rok", "Průtok") > str(dataNiger) 'data.frame': 51 obs. of 2 variables: $ Rok : int 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 ... $ Průtok: num 39.8 42.9 69.1 43.7 56.7 ... > head(dataNiger) Rok Průtok 1 1907 39.793 2 1908 42.892 3 1909 69.070 4 1910 43.653 5 1911 56.700 6 1912 45.990 Příkazem attach zpřístupníme proměnné v datovém rámci dataNiger a data vykreslíme. > attach(dataNiger) Stochastické modely časových řad (2. cvičení) 7 > plot (Rok, Průtok, main = TXT, cex.main = 0.85, type = "1") Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 ~~I-1-1-1-1— 1910 1920 1930 1940 1950 Rok Obrázek 3: Průměrně roční průtoky v řece Niger v Coulicoure (Mali). Z grafu je patrné, že závislost průtoku na čase bude složitější a proložit grafem regresní přímku by mohlo být příliš zjednodušující. Proto zkusíme daty proložit polynom vyššího stupně. Nejprve zkusíme polynom třetího stupně, kdy budeme předpokládat, že jednotlivá pozorování můžeme popsat vztahem Yí= 130+ PiXi + I32xf + /33xf + £j, kde Yi označuje průtok vody a X{ označuje rok odpovídající danému údaji. Neznámé parametry tohoto modelu nyní odhadneme. > model3 <- lm(Prutok ~ Rok + I(Rok~2) + I(Rok~3)) > summary(mode!3) Call: lmCformula = Průtok ~ Rok + I(Rok~2) + I(Rok~3)) Residuals: Min IQ Median 3Q Max -22.2179 -7.0170 -0.9496 7.6885 27.1292 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.483e+07 4.903e+06 -3.024 0.00403 ** Rok 2.302e+04 7.614e+03 3.023 0.00404 ** I(Rok~2) -1.191e+01 3.941e+00 -3.022 0.00405 ** I(Rok~3) 2.055e-03 6.800e-04 3.022 0.00406 ** Signif. codes: 0 0.001 0.01 0.05 0.1 ' 1 Residual standard error: 12.14 on 47 degrees of freedom Multiple R-squared: 0.2108, Adjusted R-squared: 0. F-statistic: 4.183 on 3 and 47 DF, p-value: 0.0105 1604 8 Stochastické modely časových řad (2. cvičení) Výsledky modelu s polynomem řádu tři zakreslíme do grafu. Odhadnutou regresní křivku zobrazíme do grafu tak, že nejprve vytvoříme vektor gridt obsahující 500 hodnot v rozmezí 1907 a 1957. Pak použijeme funkci predict, pomocí níž pro všechny prvky gridt spočítáme vyrovnané hodnoty na základě odhadnutých parametrů modelu a uložíme je do objektu Pred3. Výslednou křivku zobrazíme pomocí funkce matlines, kde jako argumenty zadáme gridt a Pred3. > n <- length(Rok) > gridt <- seq(Rok[l], Rok Lni, length.out = 500) > Pred3 <- predict(model3, data.írame(Rok = gridt)) > par (mar = c (2, 2, 1, 0) + 0.5) > plot (Rok, Průtok, type = "1", main = TXT, cex.main = 0.85) > matlines (gridt, Pred3, col = "red") Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 A \ / / A v 1910 1920 1930 1940 1950 Obrázek 4: Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) proložené polynomem třetího stupně. Zkusíme odhadnout jiný model s polynomem šestého stupně, budeme tudíž předpokládat, že jednotlivá pozorování lze popsat vztahem Yi = (30 + 0\Xi + (Í2XJ + /33Xi + 04xf + 0*,x\ + 0§x\ + Ei. Odhadneme parametry nového modelu. > model6 <- lm(Prutok ~ Rok + I(Rok~2) + I(Rok~3) + I(Rok~4) + I(Rok~5) + l(Rok~6)) > summary(modelS) Call: lmCformula = Průtok ~ Rok + I(Rok~2) + I(Rok~3) + I(Rok~4) + I(Rok~5) + I(Rok~6)) Residuals: Min IQ Median 3Q Max -22.2179 -7.0170 -0.9496 7.6885 27.1292 Stochastické modely časových řad (2. cvičení) 9 Coefficients: (3 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) -1 ,483e+07 4 ,903e+06 -3 ,024 0 ,00403 ** Rok 2 ,302e+04 7 ,614e+03 3 ,023 0 , 00404 ** I(Rok~2) -1 ,191e+01 3 ,941e+00 -3 ,022 0 ,00405 ** I(Rok~3) 2 ,055e-03 6 ,800e-04 3 ,022 0 ,00406 ** I(Rok~4) NA NA NA NA I(Rok~5) NA NA NA NA I(Rok~6) NA NA NA NA Signif. codes: 0 0.001 0.01 0.05 \ ' 0.1 x ' 1 Residual standard error: 12.14 on 47 degrees of freedom Multiple R-squared: 0.2108, Adjusted R-squared: 0.1604 F-statistic: 4.183 on 3 and 47 DF, p-value: 0.0105 Vidíme, že odhad modelu s vysokým stupněm polynomu je numericky problematický, protože se koeficienty od čtvrtého stupně výše nepodařilo spočítat. V takových případech se doporučuje proměnnou Rok centrovat a odhadnout parametry pro model s touto centrovanou proměnnou. Centrování znamená, že od každé hodnoty proměnné Rok odečteme průměr všech hodnot. Model se tak změní do podoby Yi = Po + Pí* + fozj + foz\ + p4zf + p5zf + p6zf + eh kde Zi = Xi — x. Provedeme nový odhad modelu s centrovanou proměnnou. > delta <- mean(Rok) > RokC <- Rok - delta > modeiec <- lm(Prutok ~ RokC + I(RokC~2) + I(RokC~3) + I(RokC~4) + I(RokC~5) + I(RokC~6)) > summary(modelSC) Call: lm(formula = Průtok ~ RokC + I(RokC~2) + I(RokC~3) + I(RokC~4) + I(RokC~5) + I(RokC~6)) Residuals: Min IQ Median 3Q Max -19.3877 -5.4541 -0.9673 5.1213 21.0478 Coefficients: Estimate Std. Error t value PrOltl) (Interc :ept) 6 .469e+01 2, 872e+00 22.523 < 2e-16 *** RokC -1 .757e+00 3, 909e-01 -4.494 5.02e-05 *** KRokC" 2) -2 .385e-01 5, 367e-02 -4.444 5.90e-05 *** KRokC" 3) 1 .045e-02 2, 371e-03 4.408 6.63e-05 *** KRokC" 4) 8 .809e-04 2, 268e-04 3.884 0.000341 *** KRokC" 5) -1 .165e-05 3, 209e-06 -3.631 0.000733 *** KRokC" 6) -8 .461e-07 2, 526e-07 -3.350 0.001666 ** Signif. codes: 0 x ***' 0. 001 0.01 x* ' 0.05 \ ' 0 Residual standard error: 9.36 on 44 degrees of freedom Multiple R-squared: 0.5608, Adjusted R-squared: 0.5009 F-statistic: 9.364 on 6 and 44 DF, p-value: 1.278e-06 10 Stochastické modely časových řad (2. cvičení) Výsledky opět znázorníme graficky. > gridtC <- gridt - delta > PredS <- predict(modelSC, data, frame(RokC = gridtC)) > par (mar = c (2, 2, 1, 0) + 0.5) > plot (Rok, Průtok, type = "1", main = TXT, cex.main = 0.85) > matlines(gridt, PredS, col = "red") Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 1910 1920 1930 1940 1950 Obrázek 5: Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) proložené polynomem šestého stupně. Příklad 3: Trigonometrická regrese Máme k dispozici údaje o průměrných měsíčních teplotách v New Yorku. Zjištěné teploty ve stupních Celsia za období leden 1949 - prosinec 1959 jsou uloženy v souboru nytemps. dat. Soubor obsahuje jediný sloupec se zjištěnými teplotami a není v něm uvedena hlavička. Data načteme obvyklým způsobem pomocí funkce read.table, zobrazíme strukturu vytvořeného datového rámce a vypíšeme začátek datového souboru. Proměnnou v datovém rámci pojmenujeme a zpřístupníme pomocí příkazu attach. > fileDat <- paste (data.library, > data <- read.table(fileDat) > str(data) "nytemps. dat", sep = "") 'data, frame': 168 obs. of 1 variable: $ VI: imm 11.5 11 14.4 14.4 16.5 ... > head(data) VI 11.506 11.022 14.405 14.442 16.524 17.918 Stochastické modely časových řad (2. cvičení) 11 > names(data) <- "teploty" > attach(data) Vývoj teplot v čase znázorníme graficky. > TXT <- "Průměrné mesicni teploty v~New Yorku" > n <- length(teploty) > Time <- l:n > plot(1946 + (Time - D/12, teploty, type = "1", xláb = "Cas", yláb = "Teplota") —i-1-1-1-1-1-1-1— 1946 1948 1950 1952 1954 1956 1958 1960 Cas Obrázek 6: Průměrné měsíční teploty v New Yorku v letech 1946-1959. Z grafu je patrné, že vývoj teplot v čase je periodický. Z povahy věci můžeme předpokládat, že perioda má délku 12 měsíců. Proto jako regresní funkci použijeme takovou funkci, která je periodická s periodou 12. K tomu jsou nejvhodnější funkce sinus a cosinus. Jednotlivá pozorování popíšeme rovnicí Yi = Po + /3i sin (y^xí) + cos + £i' kde Yi označuje teplotu odpovídající i-tému pozorování a X{ je odpovídající čas. K odhadu parametrů modelu opět použijeme funkci lm. > modelPeriod <- lm(teploty ~ 1 + I(sin(2 * pi/12 * Time)) + I(cos(2 * pi/12 * Time))) > summary(modelPeriod) Call: lm(formula = teploty ~ 1 + I(sin(2 * pi/12 * Time)) + I(cos(2 * pi/12 * Time))) Residuals: Min IQ Median 3Q Max -1.6248 -0.4007 0.0111 0.3264 5.4562 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 15.33289 0.05203 294.70 <2e-16 *** I(sin(2 * pi/12 * Time)) -2.53766 0.07358 -34.49 <2e-16 *** I(cos(2 * pi/12 * Time)) -2.92704 0.07358 -39.78 <2e-16 *** 12 Stochastické modely časových řad (2. cvičení) Signif. codes: 0 0.001 0.01 0.05 \ ' 0.1 x ' 1 Residual standard error: 0.6744 on 165 degrees of freedom Multiple R-squared: 0.9438, Adjusted R-squared: 0.9431 F-statistic: 1386 on 2 and 165 DF, p-value: < 2.2e-16 Nakonec zobrazíme data spolu s odhadnutou regresní krivkou. > n <- length(data$teploty) > gridt <- seq(Time[l], Time [n], length, out = 500) > PredPeriod <- predict (modelPeriod, data.frame (Time = gridt)) > par (mar = c (2, 2, 1, 0) + 0.5) > plot(1946 + (Time - D/12, data$teploty, type = "1", main = TXT, cex.main = 0.85) > matlines(1946 + (gridt - D/12, PredPeriod, col = "darkgreen") Prumerne mesicni teploty v New Yorku I-1-1-1-1-1-1-1— 1946 1948 1950 1952 1954 1956 1958 1960 Obrázek 7: Průměrné měsíční teploty v New Yorku s trigonometrickou regresní funkcí. příklad 4: jednofaktorová ANOVA Byl proveden experiment, v němž se do pastí odchytával jistý druh můry. Pasti obsahovaly vždy jeden ze tří různých typů návnady a byly umístěny v různých výškách stromů. Data z experimentu jsou uložena v souboru moths.csv a informace o datovém souboru jsou k dispozici v souboru moths_info.txt. Nejprve pomocí funkce readLines načteme soubor moths_info.txt. > íileTxt <- paste (data.library, "moths_inío. txt", sep = "") > con <- file (íileTxt) > (popis <- readLines (con)) [1] "Spruce Moth Trap" [2] "" [3] "Response: number of spruce moths found in trap after 48 hours" [4] "Factor 1: Location of trap in tree (top branches, middle branches, lower branches, ground)" [5] "Factor 2: Type of lure in trap (scent, sugar, chemical) " [6] "" [7] "" > close (con) V souboru moths. csv jsou data uložená ve třech proměnných, jejichž názvy jsou uvedeny na prvním řádku a jsou odděleny středníkem. Proto použijeme funkci read, table s argumenty sep=";" a header=TRUE. Stochastické modely časových řad (2. cvičení) 13 > fileDat <- paste (data.library, "moths. csv", sep = "") > data <- read.table(fileDat, header = TRUE, , sep = ";") > str(data) 'data.frame': 60 obs. of 3 variables: $ Location: Factor w/ 4 levels "Ground","Lower: 4444444444 ... $ Lure : Factor w/ 3 levels "Chemical","Scent",..: 2222233333 ... $ Number : int 28 19 32 15 13 35 22 33 21 17 ... Zajímá nás, jaký vliv na počet můr chycených do pasti má umístění pasti. Úrovně této proměnné typu faktor zjistíme příkazem levels. > levels(data$Location) [1] "Ground" "Lower" "Middle" "Top" Umístění pasti je v našem případě kategóriami proměnná se čtyřmi úrovněmi Ground, Lower, Middle a Top. Přibližnou představu o rozdílech mezi zjištěnými údaji pro různá umístění pastí můžeme získat z krabicového grafu. > plot(data$Location, data$Number, ylab = "Number") E CM Ground Lower Middle Top Obrázek 8: Počet chycených můr v závislosti na umístění pasti. Budeme předpokládat, že střední hodnota počtu odchycených můr se liší v závislosti na tom, v jaké výšce byla past umístěna, což můžeme zapsat E(Ygr(run([) — fl + Clground -^(^řouier) — A* ^lower E{Ymiddle) = H + (^middle E{Ytop) = /i + Oítop 14 Stochastické modely časových řad (2. cvičení) kde proměnná Y označuje počet odchycených můr a index označuje umístění pasti. Tomu odpovídá model, v němž jednotlivá pozorování popíšeme rovnicí Yjk = H + aj + ejk, (1) kde index j = 1, 2, 3,4 označuje umístění pasti odpovídající danému údaji a index k označuje pořadí pozorování mezi pozorováními se stejným umístěním pasti. Tento model se nazývá analýza rozptylu jednoduchého třídění, zkráceně jednofaktorová ANOVA. K odhadu koeficientů /x a a j použijeme opět funkci lm. > anova.1 <- lm(Number ~ Location, data.) > summary(anoval) Call: lm(formula = Number ~ Location, data = data) Residuals: Min 1Q Median 3Q Max -19.000 -4.517 -1.200 5.950 13.000 Coefficients: Estimate Std. (Intercept) 19.067 LocationLower 14.267 LocationMiddle 11.933 LocationTop 4.267 Error t value Pr(>|t|) 1.970 9.676 1.49e-13 *** 2.787 5.120 3.90e-06 *** 2.787 4.282 7.32e-05 *** 2.787 1.531 0.131 Signif. codes: 0 0.001 0.01 0.05 \ ' 0.1 Residual standard error: 7.632 on 56 degrees of freedom Multiple R-squared: 0.3779, Adjusted R-squared: 0.3446 F-statistic: 11.34 on 3 and 56 DF, p-value: 6.459e-06 Systém R ale pracuje s jinou parametrizací, než jaká je uvedena v rovnici (1). První skupina pozorování je zvolena jako referenční (zde jsou to pozorování, u nichž má proměnná Location úroveň Ground) a je zavedena následující omezující podmínka ^1 ^ground 0- Odhadnuté parametry pak mají tento význam: /x vyjadřuje střední hodnotu Y pro první (referenční) úroveň kategóriami proměnné, a j vyjadřuje rozdíl ve střední hodnotě Y mezi první a j-tou úrovní kategóriami proměnné. Z výpisu informací o modelu anoval zjistíme tyto údaje: (Intercept) = E{Yground) = 19.067 LocationLower = E(Ylower) - E(Yground) = 14.267 LocationMiddle = E(Ymiddle) - E(Yground) = 11.933 LocationTop = E(Ytop) - E(Yground) = 4.267 Stochastické modely časových řad (2. cvičení) 15 A odtud můžeme dopočítat střední hodnoty pro jednotlivé úrovně proměnné Location: E(Ygroundtk) = (Intercept) = 19.067 ĚiYiower^k) = (Intercept) + LocationLower = 19.067+ 14.267 = 33.334 Ě(Ymiddie)k) = (Intercept) + LocationMiddle = 19.067 + 11.933 = 31 Ě(Ytop :k) = (Intercept) + LocationTop = 19.067 + 4.267 = 23.334 příklad 5: dvoufaktorová ANOVA Pro stejná data týkající se odchycených můr zkusíme odhadnout model, v němž bychom zohlednili jak vliv umístění pasti, tak vliv zvolené návnady na počet můr, které se podařilo lapit do pasti. Abychom si udělali představu, jak počet chycených můr ovlivňuje návnada, vykreslíme krabicové graf pro zjištěná dat zvlášť pro jednotlivé typy návnady. > plot (data$Lure, data$Number, yláb = "Number") Chemical Scent Sugar Obrázek 9: Počet chycených můr v závislosti na typu návnady. Situaci se pokusíme popsat modelem, do něhož zahrneme vliv obou kategorických proměnných. Nejprve zkusíme odhadnout jednodušší model j = 1,2,3,4 Yjki = V + ot j + fa +£jkl, kde k = 1,2,3 l — 1,..., n j k kde aj zachycuje vliv umístění pasti a zachycuje vliv fc-tého typu návnady a pro počty pozorování platí Yľj=i Sfc=i njk = n = 60. > anova2 <- lm(Number ~ Lure + Location, data) > summary (anova2) 16 Stochastické modely časových řad (2. cvičení) Call: lm(formula = Number ~ Lure + Location, data = data) Residuals: Min 1Q Median 3Q Max -17.4500 -5.1208 -0.5167 6 0625 12.9333 Coefficients: Estimate Std. Error t value PrOltl) (Intercept) 19.883 2 415 8.234 4.14e-ll *** LureScent -2.750 2 415 ■ -1.139 0.260 LureSugar 0.300 2 415 0.124 0.902 LocationLower 14.267 2 788 5.117 4.23e-06 *** LocationMiddle 11.933 2 788 4.280 7.70e-05 *** LocationTop 4.267 2 788 1.530 0.132 Signif. codes: 0 "***' 0.001 ' **' 0 01 x* 0.05 0. Residual standard error: 7.636 on 54 degrees of freedom Multiple R-squared: 0.3995, Adjusted R-squared: 0.3439 F-statistic: 7.184 on 5 and 54 DF, p-value: 3.236e-05 Systém R opět zvolil první úrovně faktorů jako referenční. Proto koeficient označený v tabulce jako (Intercept) představuje odhadnutou střední hodnotu chycených můr pro Location = Ground a Lure = Chemical. Zbylé koeficienty představují odchylky středních hodnot od referenční v případě, že se změní úroveň jedné z proměnných. V modelu anova2 jsme mlčky předpokládali, že vliv proměnných Location a Lure se sčítá, nebrali jsme v úvahu možnost, že by mohlo docházet ke složitějším interakcím mezi nimi. Jak se mění průměrné hodnoty chycených můr pro různé kombinace úrovní obou vysvětlujících proměnných si můžeme prohlédnout v následujícím grafu. > interaction.plot(data$Location, data$Lure, data$Number, xlab = "Location", ylab = "Number", trace.label = "Lure") Obrázek 10: Porovnání průměrného počtu chycených můr v závislosti na umístění pasti a typu návnady. Stochastické modely časových řad (2. cvičení) 17 Nyní odhadneme složitější model, který bude mít tvar j = 1,2,3,4 Yjki = V + otj + Pk +(otP)jk + £jki, kde k = 1,2,3 l — 1, . . . , Tljk kde a j zachycuje vliv umístění pasti, /% zachycuje vliv fc-tého typu návnady a {a(i)jk zachycuje vliv interakce mezi j-tou úrovní proměnné Location a fc-tou úrovní proměnné Lure. > anova.3 <- lm(Number ~ Lure * Location, data) > summary(anova3) Call: lm(formula = Number ~ Lure * Location, data = data) Residuals: Min 1Q Median 3Q Max -15.80 -5.00 -0.90 5.85 14.20 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 19. 200 3. 555 5. 400 2.04e-06 *** LureScent -2. 200 5. 028 -0. 438 0. 66367 LureSugar 1. 800 5. 028 0. 358 0. 72191 LocationLower 16. 800 5. 028 3. 341 0. 00162 ** LocationMiddle 12. 600 5. 028 2. 506 0. 01565 * LocationTop 3. 800 5. 028 0. 756 0. 45347 LureScent:LocationLower -1. 000 7. 111 -0. 141 0. 88875 LureSugar:LocationLower -6. 600 7. 111 -0. 928 0. 35795 LureScent:LocationMiddle -1. 800 7. 111 -0. 253 0. 80124 LureSugar:LocationMiddle -0. 200 7. 111 -0. 028 0. 97768 LureScent:LocationTop 0. 600 7. 111 0. 084 0. 93310 LureSugar:LocationTop 0. 800 7. 111 0. 113 0. 91089 Signif. codes: 0 0. 001 v **' 0.01 ' *' 0. 05 0.1 x ' 1 Residual standard error: 7.95 on 48 degrees of freedom Multiple R-squared: 0.4214, Adjusted R-squared: 0.2888 F-statistic: 3.178 on 11 and 48 DF, p-value: 0.002653 příklad 6: ANCOVA Při experimentu byla zkoumána závislost rychlosti růstu populace studenokrevných organismů na teplotě prostředí. V laboratoři byla sledována rychlost populačního růstu (proměnná rate) u dvou druhů roztočů (kategoriální proměnná genus se dvěma úrovněmi genA a genB) při různých teplotách (proměnná temp). Data jsou obsahem souboru mite.txt, v němž je na prvním řádku uvedena hlavička a hodnoty pro jednotlivé proměnné jsou odděleny tabelátorem. Proto k načtení použijeme funkci read.table s argumentem sep = "\t". > íileDat <- paste (data, library, "mite, tjct", sep = "") > roztoči <- read, table (ÍileDat, sep = "\t", header = TRUE) > str(roztoči) 18 Stochastické modely časových řad (2. cvičení) 'data.frame': 22 obs. of 3 variables: $ temp : num 10 12.5 15 18 20 22.5 25 27.5 30 32.5 ... $ rate : num 0.01 0.0759 0.1083 0.1505 0.1774 ... $ genus: Factor w/ 2 levels "genA","genB": 1111111111... > attach(roztoči) Načtená data znázorníme graficky. Abychom si udělali představu o tom, jak se výsledky liší pro různé druhy roztočů, znázorníme jinou barvou údaje získané pro roztoče druhu A (bílá barva) a pro roztoče druhu B (černá barva). Toho dosáhneme tak, že nejprve nakreslíme „prázdný" graf a pak do něj zvlášť vykreslíme body reprezentující data pro druh A a pro druh B. > plot (temp, rate, type = "n") > points (temp [genus == "genA"], rate [genus == "genA"]) > points (temp [genus == "genB"], rate [genus == "genB"], pch = 16) > legend("topieft", legend = cC'genA", "genB"), pch = c(l, 16), bty = "n") temp Obrázek 11: Závislost rychlosti růstu na teplotě pro dva druhy roztočů: "genA" (bíle), "genB" (černě). Podle obrázku bychom mohli usoudit, že závislost rychlosti růstu populace na teplotě není lineární a navíc pro druh A je většinou o něco vyšší. To nás vede k tomu, že model popisující vztah proměnných, by mohl mít následující podobu Yjk = H + aj + /3xjk + -fx2jk + ejk, kde Y označuje rychlost růstu populace, x je teplota, j je index označující druh roztoče a k označuje konkrétní pozorování v rámci daného druhu. Modely, v nichž figuruje alespoň jedna spojitá a alespoň jedna kategóriami proměnná, se nazývají analýza kovariance (ANCOVA). > ancova <- ImCrate ~ temp + I(temp'2) + genus) > summary(ancova) Stochastické modely časových řad (2. cvičení) 19 Call: lm(formula = rate ~ temp + I(temp~2) + genus) Residuals: Min 1Q Median 3Q Max -0.013530 -0.006213 -0.002241 0.003962 0.017629 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.111e-01 1.779e-02 -17.489 9.63e-13 *** temp 4.060e-02 1.694e-03 23.966 4.15e-15 *** I(temp~2) -8.154e-04 3.719e-05 -21.927 1.96e-14 *** genusgenB -1.224e-02 4.128e-03 -2.966 0.00828 ** Signif. codes: 0 "***' 0.001 "**' 0.01 "*' 0.05 ".' 0.1 x ' 1 Residual standard error: 0.009682 on 18 degrees of freedom Multiple R-squared: 0.9753, Adjusted R-squared: 0.9712 F-statistic: 237.1 on 3 and 18 DF, p-value: 1.186e-14 Podobně jako v předchozích příkladech byl druh A vybrán jako referenční, a proto odhadnutý parametr genusgenB vyjadřuje odchylku druhu B od druhu A. Výsledky odhadu modelu zobrazíme do grafu. > grid <- seq(from = min(temp), to = max(temp), length.out = 300) > predA <- predict (ancova, data, frame (temp = grid, genus = repC'genA", 300))) > predB <- predict (ancova, data, frame (temp = grid, genus = repC'genB", 300))) > plot (temp, rate, type = "n") > points (temp [genus == "genA"] , rate [genus == "genA"]) > points (temp [genus == "genB"], rate [genus == "genB"], pch = 16) > matlines(grid, predA, col = "navy") > matlines(grid, predB, col = "darkgreen") > legend(x = 10.5, y = 0.2, legend = cC'genA", "genB"), pch = c(l, 16), col = c("navy", "darkgreen"), bty = "n") > legend(x = 9, y = 0.18, legend = cC'genA", "genB"), lty = 1, col = c("navy", "darkgreen"), bty = "n") 20 Stochastické modely časových řad (2. cvičení)