RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 1 M5201 - 2. CVIČENÍ: Regresní analýza v R 1 Uvod 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é (kovariáty, regresory) 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é Y1,Y2,...,Yn a jim odpovídající hodnoty vysvětlujících proměnných Xll, X12, • • • , Xin, X21, X22i • • • 1 X2rn • • • > Xmi, . . . , Xmn a pomocí regresního modelu můžeme popsat jednotlivá pozorování Yi — f(x{±, X{2) • • • ) X{m, £j), i — l,...,fl Nejjednodušším případem regresního modelu je klasický 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ětlujících proměnných s neznámými regresními koeficienty. Pomocí rovnice můžeme lineární regresní model zapsat takto Yi = Po + Plfl(xn,Xi2, . . . ,Xim) + P2Í2(xn,Xi2, . . .,Xim) -\-----h Pkfk (xil, Xi2, . . . , Xim) + £j, kde pro i = 1, Y Xil) • • • ) Xim fl, f2, ■ ■ ■ , fk ßo,ßi, ■ ■ ■ ,ßk sou vysvětlované proměnné sou vysvětlující proměnné (kovariáty, regresory) sou známé (předem zvolené) funkce vysvětlujících proměnných sou neznámé regresní koeficienty £j jsou nekorelované náhodné veličiny s nulovou střední hodnotou a konstat- ní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. Stochastické modely časových řad (2. cvičení) 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í ~ IID(0,a2). Pro testování hypotéz je třeba přidat předpoklad, že rozdělení chybové složky je normální. Pro odhad neznámých parametrů regresních modelů existuje více metod. Pro lineární regresní modely se obvykle používá metoda nejmenštch čtverců. Odhady parametrů se obvykle označují stříškou, např. Po,Pi,.... Hodnoty vysvětlované proměnné Yi, které získáme po dosazení odhadnutých parametrů do regresní rovnice Yi = Po + Plfl{xn,Xi2, Xim) + P2Í2{xil,Xi2, ■■■ , Xim) H-----h Pkfk{xil,Xi2, ... , Xim) 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ě. Poznámka Při maticovém popisu regresního modelu obvykle provedeme určité přeznačení systematická složka Y=p0+Pifi( Zil> %i2) • • • ) %ip ) +@2 Í2{zn,Zi2, . . . , Zip) + ■ ■ ■ + Pm fm(Zil, zí2, ■ ■ ■ , Zim) +£'í označíme xn označíme Xi2 označíme Xim takže pracujeme s modelem Yl = Po + PlXll H-----VPmXik + Ei (Yi tj. maticově Yn — Po + P\Xn\ + • • • + PmXrim ~\~ &n \Yri Y X(matice plánu) Uvažujme klasický regresní model plné hodnosti, kde chyba má normální rozdělení: Y = X/3 + e A h(X) = h(X'X) =m + l A n > m + 1 A e~ Nn(0, a2In) Odhad neznámých parametrů /3 provedený metodou nejmenštch čtverců je řešením normálních rovnic X'X/3 = X'Y a platí 3 = (X'X)_1X'Y. Označme Ý = X/3 = X (X'X)-1 X'Y = HY H e = Y - Ý = (I - H)Y = MY = M(X/3 + e)) = MX/3 + Me = (I - H)e M =0 s2 = ^T = ^^(Y-Ý)/(Y-Ý) = ^^é/e = ^^Y(I-H)Ý = ^^e/(I-H)e n—p—1 n—m—l^ > v > n—m—l n—m—l v > n—m—l v > RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova, Pak připomeňme následující vlastnosti klasického regresního modelu (plné hodnosti) ■ E(3 = (3 (nestrannost odhadu) ■ Es2 = nE^-i = °'2' tj. Je nestranným odhadem rozptylu . Y ~ Nn(Xp,a2In) m e~7Vn(0,a2(I-H)) . 3~iVm+1(/9,<72(X'X)-1) ■ ~f ~ X2(n — m — 1) ■ /3 a s2 jsou stochasticky nezávislé Tj = AJL ~ í(n-m-l), kde (X'X)"1 = (%)m=0,. ^ = ^(/32-/32)/W-1(32-/32) ~ F(q,n-m-l), kde (X'X)-1=(^), ^(l;) a fc(W) = í T= /2c^~tw ~ŕ(n-m-l), kde c = (c0, ci,..., cm)' a E(dp)=c'p ^ y.-y, ~iv(o,a2(i + x (x'x)-1^)) kde x,' = (xjo, • • • *Xim) Je ^-tý řádek matice plánu X Na základě předchozích vlastností lze odvodit následující intervaly spolehlivosti 4 Stochastické modely časových řad (2. cvičení) Intervaly spolehlivosti pro parametry (5j 3 =0, ...,p (Pj - *i-f {n-p-V)Sy/vj], j3j +íi-f (n-p-l)s^u~) pro střední hodnotu predikce £rí = £x-/3 = x-/3 x^3-Í!_f (n-m-l)s1/x^(X'X)-1xi, x^/3 + íi_f (n-m-l)s1/x^(X'X)-1xi pro predikci £ = x$ i = 1,..., n x^-t^n-m-ljST/l+x^X'Xj-ix^x^í+íi-f (n-m-l)s1/l+x^(X'X)-1xi kde í1_a.(n —m —1) je 1 — kvantil Studentova rozdělení o n —m —1 stupních volnosti 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 Y{ = [3q + f3\Xn + $2X21 + £«, zatímco y ~ I(xl + x2) definuje model Yi = (3q + f3i(xu + X2i) + £j. Konkrétní příklady zadávání modelů v systému R jsou ukázány v následující tabulce. RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 5 Formální zápis Syntaxe v R Popis Yi = n + el y ~ i Model pouze s absolutním členem Yi = ßo + ßi%i + £i 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 \og{Yi) = ßo + ßixi + el log (y) - x Lineární model se spojitou vysvětlující x a s logaritmicky transformovanou vysvětlovanou proměnnou Yi = ßo + ßlXi + ßix\ + Ei y ~ y ~ x + I(x' poly(x, •2), 2) Kvadratický model se spojitou vysvětlující proměnnou x (regresní parabola) Yi= ß0+ ßlXil + ßlXa + £i 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 = ßo + ßxXiX + ß2xl2 + ß3XtlXi2 + et y ~ y ~ xl*x2, xl+x2+I(xl*x2) Model se dvěma spojitými vysvětlujícími proměnnými xi & x2 s křížovým členem Yij — /i + Ui + Sij y ~ A Model s jedním faktorem (jed-nofaktorová analýza rozptylu, ANOVA) Yijk = n + a>i + ßj + Eijk y ~ A + B Model se dvěma faktory bez interakce (dvoufaktorová analýza rozptylu) Yijk = n + oti + ßj + {aß) i j + eljk y ~ A + B + A:B, y ~ A * B Model se dvěma faktory s interakcí (dvoufaktorová analýza rozptylu) Yij m — ^ + ai + ßj + 7fc + EijM y ~ A + B + C Model se třemi faktory bez interakcí (třífaktorová analýza rozptylu) Yij m = M + ®i + ßj + 1k + {aß)ij + («7)ifc + {ßl)jk + {aßrfijk + Sij m y ~ y ~ A+B+C+A: A*B*C :B+A: :C+B:C+A:B:C Model se třemi faktory s interakcemi (třífaktorová analýza rozptylu) Yijki = M + cti + ßj + 7fc + {aß) i j + {a^)ik + {ßl)jk + £ijki y ~ y ~ A+B+C+A: (A+B+C)' :B+A: "2 :C+B:C Model se třemi faktory obsahující kromě hlavních efektů maximálně dvojné interakce (třífaktorová analýza rozptylu) - j-Jf 1 J 1 ß X y£J 1 SXJ y ~ x+A Model se spojitou proměnnou x a faktorem A bez interakce interakce (ANCOVA) — H~ Oĺj ~h ßXij ~\~ öjXij ~h Sij 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ů (spojité veličiny značíme písmeny z konce abecedy: x, Y,..., faktory písmeny ze začátku abecedy: A,B,...) 6 Stochastické modely časových řad (2. cvičení) 3 Příklady Příklad 1: Regresní přímka Máme k dispozici údaje o počtu australských domácností připojených v letech 1998-2008 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) O O o CD O O O m o o o o o ro E o "O o _)_i o CD o O co O Q. o o o OJ o o o "1-ľ 1998 2000 2006 2008 2002 2004 Rok Obrázek 1: Počet domácností s připojením na internet v letech 1998-2008 Budeme uvažovat regresní model, ve kterém počet domácností s připojením na internet závisí lineárně na čase. Model můžeme popsat regresní rovnicí Yi = f3o + (3\Xi + eí 1, - - -, 11 kde RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 7 Yi 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)> /?i 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 summaryO. > modeli <- Im (domácnosti ~ time) > summary (modeli) Call: lm(formula = domácnosti ~ time) Residuals: Min IQ 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 ' ' 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 p0 = -948407.45 (Intercept) ft = 475.36 (time) 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. Do grafu navíc dokreslíme interval spolehlivosti kolem střední hodnoty a predikční interval spolehlivosti (který je širší a měl by obsahovat cca 95% dat). > n <- length(time) > gridt <- seq (time [ 1 ], time fuj, length, out = 500) > ci.coní <- predict(modell, newdata = data.frame(time = gridt), interval = "confidence") > ci.pred <- predict (modell, newdata = data, frame (time = gridt), interval = "prediction") > yrange <- range (c (domácnosti, ci.predl, 2], ci.predl, 3D) > plot (time Lc(l, n)l, yrange, type = "n", xlab = TxtX, ylab = TxtY) > matlines (gridt, ci.confL, 1], lty = 1, lwd = 2, col = "red2") > matlines (gridt, ci.confL, 2:3], lty = 2, lwd = 1.5, col = "dodgerblue") > matlines (gridt, ci.predL, 2:3], lty = 3, lwd = 1.5, col = "darkgreen") > points(time, domácnosti, pch = 21, cex = 0.85, bg = "redl") 8 Stochastické modely časových řad (2. cvičení) o o o o o 1998 2000 2002 2004 2006 2008 Rok Obrázek 2: Grafické vyjádření regresní přímky pro data Počet domácností s připojením na internet v letech 1998-2008 spolu s intervaly spolehlivosti kolem střední hodnoty a s predikčním intervalem Pro zajímavost totéž provedeme v trochu jiné grafické úpravě. > plot (time [c (1, n)], yrange, type = "n", xlab = TxtX, ylab = TxtY) > xx <- cCgridt, rev(gridt)) > yy <- c(ci . coní [, 2], rev(ci . coní [, 3])) > polygon(xx, yy, col = "gray75", border = "gray75") > yy <- c(ci. coní[, 2], rev(ci .pred[, 2])) > polygon(xx, yy, col = "gray85", border = "gray85") > yy <- c (ci . coní[, 3], rev(ci .předl, 3D) > polygon(xx, yy, col = "gray85", border = "gray85") > matlines(gridt, ci.coníL, 11, Ity = 1, lwd = 2, col = "red2") > points (time, domácnosti, pch = 21, cex = 0.85, bg = "yellowgreen") RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 9 ^-1-1-1-1- 1998 2000 2002 2004 2006 2008 Rok Obrázek 3: Grafické vyjádření regresní přímky pro data Počet domácností s připojením na internet v letech 1998-2008 spolu s intervaly spolehlivosti kolem střední hodnoty a s predikčním intervalem - s využitím ploch 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. f rame': 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 ... 10 Stochastické modely časových řad (2. cvičení) > 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) > 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-r~ 1910 1920 1930 1940 1950 Rok Obrázek 4: 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 Yi = po + (3ixi + fox? + /33xf + eí, 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. RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 11 > model3 <- lm(Prutok ~ Rok + I(Rok~2) + I(Rok~3)) > summary (model3) Call: lm(formula = Průtok ~ Rok + I(Rok~2) + I(Rok~3)) Residuals: Min 1Q 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.1604 F-statistic: 4.183 on 3 and 47 DF, p-value: 0.0105 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 na základě odhadnutých parametrů modelu vyrovnané hodnoty jednak odhadnuté hodnoty, kromě toho také intervaly spolehlivosti kolem střední hodnoty a predikční intervaly > n <- length(Rok) > gridt <- seq(Rok[l], Rok[n], length.out = 500) > ci.conf <- predict(model3, newdata = data.frame(Rok = gridt), interval = "confidence") > ci.pred <- predict (model3, newdata = data, frame (Rok = gridt), interval = "prediction") > yrange <- range(c(Prutok, ci.pred[, 2], ci.pred[, 3])) > par (mar = c(2, 2, 1, 0) + 0.5) > plot(Rok[c(l, n)l, yrange, type = "n", main = TXT, cex.main = 0.85) > matlines(gridt, ci.conf[, 11, lty = 1, lwd = 2, col = "red2") > matlines (gridt, ci.conf [, 2:3], lty = 2, lwd = 1.5, col = "dodgerblue") > matlines (gridt, ci.predl, 2:3], lty = 3, lwd = 1.5, col = "darkgreen") > points(Rok, Prutok, type = "o", pch = 21, cex = 0.85, bg = "redl") 12 Stochastické modely časových řad (2. cvičení) Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 o . o — o o 00 o CD O OJ 1910 1920 1930 1940 1950 Obrázek 5: Grafické znázornění polynomické regrese 3. stupně pro data Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) spolu s intervaly spolehlivosti Opět intervaly spolehlivost zdůrazníme pomocí ploch. > par (mar = c(2, 2, 1, 0) + 0.5) > plot(Rok[c(l, n)], yrange, type = "n", main = TXT, cex.main = 0.85) > xx <- cCgridt, rev(gridt)) > yy <- c(ci . coní [, 2], rev(ci . coní [, 3])) > polygon(xx, yy, col = "gray75", border = "gray75") > yy <- c(ci. coní[, 2], rev(ci .pred[, 2])) > polygon(xx, yy, col = "gray85", border = "gray85") > yy <- c (ci . coní[, 3], rev(ci .předl, 3D) > polygon(xx, yy, col = "gray85", border = "gray85") > matlines(gridt, ci.coníL, 11, lty = 1, lwd = 2, col = "red2") > points(Rok, Průtok, type = "o", pch = 21, cex = 0.85, bg = "yellowgreen") RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 13 Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 o o o 00 o CD O o OJ 1910 1920 1930 1940 1950 Obrázek 6: Grafické znázornění polynomické regrese 3. stupně pro data Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) spolu s intervaly spolehlivosti - pomocí ploch Zkusíme odhadnout jiný model s polynomem šestého stupně, budeme tudíž předpokládat, že jednotlivá pozorování lze popsat vztahem Y = po + f3ixi + p2x\ + fox\ + p4xj + phx\ + p%x\ + £,t. Odhadneme parametry nového modelu. > model6 <- lm(Prutok ~ Rok + I(Rok~2) + I(Rok~3) + I(Rok~4) + I(Rok~5) + I(Rok~6)) > summary(modelS) Call: lm(formula = 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 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 14 Stochastické modely časových řad (2. cvičení) I(Rok~6) NA NA NA NA 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.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=/3o + fa* + /32zf + p3zf + p4zf + p5zf + pezf + eí, kde Zi = Xi — x. Provedeme nový odhad modelu s centrovanou proměnnou. > delta <- mean(Rok) > RokC <- Rok - delta > modelSC <- 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: (Intercept) RokC I(RokC~2) I(RokC~3) I(RokC~4) I(RokC~5) I(RokC~6) Estimate 6.469e+01 -1.757e+00 -2.385e-01 1.045e-02 8.809e-04 -1.165e-05 -8.461e-07 Std. Error t value 2.872e+00 22.523 3.909e-01 -4.494 5.367e-02 -4.444 2.371e-03 4.408 2.268e-04 3.884 3.209e-06 -3.631 2.526e-07 -3.350 PrOltl) < 2e-16 *** 5.02e-05 *** 5.90e-05 *** 6.63e-05 *** 0.000341 *** 0.000733 *** 0.001666 ** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 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 Výsledky opět znázorníme graficky. > n <- length(RokC) > gridt <- seq(Rok[l] , Rok Lni, length.out = 500) > ci.coní <- predict(modelSC, newdata = data.frame(RokC = gridt - delta), interval = "confidence") RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 15 > ci.pred <- predict (modelSC, newdata = data.frame(RokC = gridt - deltas, interval = "prediction "J > yrange <- range(c(Prutok, ci.pred[, 2], ci.pred[, 3])) > par (mar = c(2, 2, 1, 0) + 0.5) > plot(Rok[c(l, n)], yrange, type = "n", main = TXT, cex.main = 0.85) > matlines (gridt, ci.confL, 1], lty = 1, lwd = 2, col = "red2") > matlines (gridt, ci.confL, 2:3], lty = 2, lwd = 1.5, col = "dodgerblue") > matlines (gridt, ci.predi, 2:3], lty = 3, lwd = 1.5, col = "darkgreen") > points(Rok, Prutok, type = "o", pch = 21, cex = 0.85, bg = "redl") Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 o 00 o CD O o OJ 1910 1920 1930 1940 1950 Obrázek 7: Grafické znázornění polynomické regrese 6. stupně pro data Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) spolu s intervaly spolehlivosti Opět intervaly spolehlivost zdůrazníme pomocí ploch. > parGnar = c(2, 2, 1, 0) + 0.5) > plot(Rok[c(l, n)], yrange, type = "n", main = TXT, cex.main = 0.85) > xx <- c(gridt, rev(gridt)) > yy <- c(ci . conf [, 2], rev(ci . conf [, 3])) > polygon(xx, yy, col = "gray75", border = "gray75") > yy <- c(ci. conf [, 2], rev(ci .pred[, 2])) > polygon(xx, yy, col = "gray85", border = "gray85") > yy <- c (ci . conf [, 3], rev(ci .předl, 3])) > polygon(xx, yy, col = "gray85", border = "gray85") > matlines (gridt, ci.confL, 1], lty = 1, lwd = 2, col = "red2") > points(Rok, Prutok, type = "o", pch = 21, cex = 0.85, bg = "yellowgreen") 16 Stochastické modely časových řad (2. cvičení) Průměrné rocni průtoky vody v rece Nigeru v Coulicoure (Mali) v letech 1907 az 1957 o o OJ 1910 1920 1930 1940 1950 Obrázek 8: Grafické znázornění polynomické regrese 6. stupně pro data Průměrné roční průtoky vody v řece Niger v Coulicoure (Mali) spolu s intervaly spolehlivosti - pomocí ploch Na závěr příkladu nesníme zapomenout zrušit zpřístupnění jména proměnných v datovém rámci > detach(dataNiger) 4 Úkoly 1. Načtěte datový soubor japanese-vehicle. txt obsahující údaje o množství motorových vozidel vyrobených v Japonsku v letech 1947-1989. - Data vykreslete do grafu. - Odhadněte parametry následujících čtyř modelů: Model 1 Model 2 Model 3 Model 4 Yi = Po + fiiXi + Si Yí = Po + PiXi + fox i + Si Yi = Po + PlXi + P2XJ + p3X^ + Ei Yi = Po + PlXi + p2xj + p3X^ + p4xf + Ei, kde Yi je počet vyrobených vozidel a X{ je čas. Pokud to bude vhodné, pracujte s centrovaným časem. Pro každý model do grafu vykreslete data, odhadnutou regresní křivku, interval spolehlivosti kolem střední hodnoty a predikční interval spolehlivosti. RNDr. Marie Forbelská, Ph.D., Mgr. Marie Levákova 17 - Na závěr vykreslete grafy pro všechny modely do jednoho obrázku. 2. Do datového souboru US_population.txt byly zaznamenány údaje o počtu obyvatel v USA v desetiletých intervalech v rozmezí let 1790-1980. - Data načtěte. - Do grafu vykreslete (a) původní hodnoty časové řady (b) časovou řadu, v níž jsou počty obyvatel zlogaritmovány. - Odhadněte parametry v modelech: kde Yi je počet obyvatel a xí je rok. - Pro každý model do grafu vykreslete data, odhadnutou regresní křivku, interval spolehlivosti kolem střední hodnoty a predikční interval spolehlivosti. U modelu s logaritmem vykreslete jak graf se zlogaritmovanými hodnotami, tak s hodnotami na původní škále. - Na závěr vykreslete grafy pro oba modely do jednoho obrázku (v původním měřítku bez logaritmování). Model 2: Model 1: Yi = ßo + ßiXi + ß^x\ + Ei In Yi = ßo + ß\Xi + ß2x2 + Ei,