M5120 - cvičení Lineární regresní model _i Ondřej Pokora (pokora@math.muni.cz) Ustav matematiky a statistiky, Přírodovědecká fakulta, Masarykova univerzita, Brno (podzim 2015) Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 1/24 Lineární regresní model (LRM) ► popisuje lineární závislost V ► Y je závislá proměnná, náhodná veličina ► X\,...,Xp jsou nenáhodné vysvětlující proměnné, tzv. regresory ► e je náhodná chyba, E(e) = 0, s neznámým konstatním rozptylem D(e) = cr2 ► /3o,j6i,.. .,/3p jsou neznámé parametry ► úkol regresní analýzy: na základě opakovaných měření závislé proměnné za různých hodnot regresorů optimálně určit parametry ]6o,j6i,... ,/3p modelu ► předpokládáme, že měření je opakováno n krát, tzn. pro i = 1, máme: yi = j8o + ^j8;^; + ei ;=i ► náhodné chyby £i,... ,£n jsou nezávislé stejně rozdělené náhodné veličiny (i.i.d.) M5120 - cvičení - Lineární regresní model 2/24 m LRM (y,Xj6,(72í) plné hodnosti v ► Yj = j6o + Y_jPjxi,j + £i v maticovém tva ru \1 "NT" X + ► V je náhodný vektor n pozorování ► regresory tvoří nenáhodnou n x k matici plánu (design matrix) X ► j6 je vektor k = p + 1 neznámých parametrů ► závislost je lineární vzhledem k parametrům j6y ► vektor náhodných chyb má kovariační matici D(e) = cr2 Jn ► n > k = p + 1, tj. počet pozorování je větší než počet parametrů ► matice plánu je plné hodnosti: h(X) = k = p + 1 Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 3/24 Metoda nejmenších čtverců (MNČ) ► V LRM Y = X/S + e chceme najít vektor parametrů /S tak, aby naměřené hodnoty Y byly optimálně aproximovány vektorem X/S. ► Metoda nejmenších čtverců (ordinary least-square method) stanovuje odhad /S jako bod minima penalizační funkce, která je součtem čtverců odchylek: ní V \2 E y* - čo - E Hi h = (y - - x/s) —► min =^ i6 i=l \ j=l ► Při splnění podmínek pro LRM plné hodnosti existuje vždy právě jedno řešení této minimalizační úlohy. To lze nalézt vyřešením soustavy normálních rovnic X7Xj6 = X'Y ► Odhad vektoru parametrů v LRM metodou nejmenších čtvreců je tedy tvaru J8 = (X'X)"1*^ Poznámka: při numerických výpočtech se inverzní matice nepočítá přímo, ale využívá se např. Q-R rozkladu. Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 4/24 Veličiny v LRM (1) ► Aproximované hodnoty závislé proměnné (fitted values, Y-hat) Y = Xj6 = X(X/X)"1X/y ► Rezidua (residuals): r = Y — Y ► Reziduálni součet čtverců z=l z=l kvantifikuje velikost variability, kterou se nepodařilo LRM vysvětlit. ► Odhad rozptylu cr2 náhodných chyb C = S - - — - n — k n — p — l ► Standardní reziduálni chyba (residual Standard error): a = s = a/Š^ Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model Veličiny v LRM (2) ► Regresní součet součet čtverců n Sr = £(Ýi-Yf i=l kvantifikuje velikost variability, kterou se LRM podařilo zachytit. Je dán součtem kvadrátů odchylek aproximovaných hodnot od výběrového průměru ► Celkový součet součet čtverců je násobkem výběrového rozptylu: n St = ^{Yi-Y)1 = {n-l)S\ i=i ► Platí: St = Sr + S, Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 6/24 Koeficient determinace v LRM ► Koeficient determinace (coefficient of determination, R-squared) ^2 ^ Se Sr Sr/n ML-odhad vysvětleného rozptylu St St St I n ML-odhad celkového rozptylu se používá ke kvantifikaci poměrné části variability, kterou se LRM podařilo vysvětlit ► R2 G [0,1] ► R2 = koeficient mnohonásobné korelace = korelační poměr ► Korigovaný koeficient determinace (adjusted R-squared, R-bar-squared) Ř2 = l-(l-R2)^-l = n — k Se/(n — k) i^i odhad rozptylu chyb St/ (n — 1) S y výběrový celkový rozptyl M5120 - cvičení - Lineární regresní model 7/24 m Řešení LRM pomocí R V R řeší LRM příkaz lm (linear model): model <- lm (formule, data=tabulka) Máme-li zvlášť vektor regresorů (x) a pozorování (V), datovou tabulku (data frame) vytvoříme příkazem: tabulka <- data.frame (x, Y) Zápis tzv. formule pro některé regresní funkce: Y = j80 + j8ix Y - - x, nebo Y - 1 + x, abolutní člen je vkládán implicitně ^^^^^^^^^ Y - - 0 + x y = fa+fa*+fax2 Y - - x + I(x~2) ^^^^^^^^^ Y - - 0 + I(abs(x)) Y = fip + fre* Y - - I(exp(x)) Y = fa + falnx Y - - Klog(x)) Detailní výsledky a další číselné charakteristiky získáme příkazem výsledky <- summary (model), příp. s parametrem correlation=TRUE pro výpočet výběrové korelační matice parametru. M5120 - cvičení - Lineární regresní model s/24 m Veličiny v LRM pomoci R MNC-odhady parametru model$coefficients coef(model) odhady, směrodatné odchylky, testy významnosti, p-hod noty vysledky$coefficients coef(vysledky) Y aproximované hodnoty model$f itted.values fitted.values(model) r rezidua model$residuals residuals(model) n —k stupně volnosti model$df.residual matice plánu model.matrix(model) s odhad směrod. odchylky chyb vysledky$sigma koeficient determinace vysledky$r.squared Ŕ2 korigovaný koef. determinace vysledky$adj.r.squared (F, k — l,n — k) celkový F-test vysledky$fstatistic n — k, k) stupně volnosti vysledky$df korelační matice pro jS vysledky$correlation Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 9/24 Grafy regresních funkcí v R Připravíme souřadnicový systém vhodných rozměrů (první dva vektory min-max), do nějž zatím nic kreslit nebudeme (type="n"), přidáme popisky: plot (c(0,70), c(-5,60), type="n", xlab="osa x", ylab="osa y") Pomocí bodů vykreslíme data z tabulky, tzn. (x, Y). První dva parametry jsou vektory x-ových a y-ových souřadnic, následují grafické parametry: points (tabulka$x, tabulka$Y, col=4, pch=24, lwd=1.5, cex=1.0) Zvolíme si dostatečně hustou síť x-ových souřadnic (x*): xx <- seq (0, 70, by=0.1) Dopočítáme k nim odpovídající y-ové souřadnice, tzn. Y*: YY <- predict (model, data.frame (x=xx)) Vykreslíme graf funkce jako křivku (x*,Y*), podobně jako body: lineš (xx, YY, col=2, lwd=1.5, lty=2) Obrázek můžeme uložit mj. příkazem dev. copy2pdf (f ile=" obrazek.pdf") M5120 - cvičení - Lineární regresní model o CD O LO O CD CD N O o > o o ^1- o co o CM A A ' o H 0 10 20 30 40 50 60 70 matka Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 11/24 LRM - příklady (1) Datové soubory k následujícím příkladům jsou dostupné na serveru v adresáři /Vyuka/R/M5120/data, odkud si je zkopírujte. Pro jednotlivé soubory řešte lineární regresní modely, spočítejte odhady parametrů a další číselné charakteristiky, a modely porovnejte. Vykreslete do jednoho obrázku data i grafy spočítaných regresních funkcí. Později přidejte testování významnosti modelu (F-test) a testy významnosti parametrů (ŕ-testy). Příklad 1 (Datový soubor C3H603.txt) Pomocí regresní přímky (a regresní paraboly) spočítejte MNČ-odhady vektoru parametrů /S, aproximace Y, reziduálni součty čtverců Se a s2 v LRM (Y, X /S, o21). X 40 64 34 15 57 45 Y 33 46 23 12 56 40 Tabulka uvádí výsledky měření závislosti množství kyseliny mléčné ve 100 ml krve u matky-provorodičky x, a jejího novorozence, Y. Obě veličiny jsou uváděny jako hmotnost v mg. M5120 - cvičení - Lineární regresní model 12/24 LRM - příklady (2) Příklad 2 (Datový soubor roztaznost.txt) Pomocí regresní přímky procházející počátkem spočítejte MNC-odhady parametrů /S, aproximace Y, reziduálni součty čtverců Se a s2 v LRM (Y,X [í,cr21) pro data X 10 20 30 40 50 60 Y 0.18 0.35 0.48 0.65 0.84 0.97 Jedná se o měření koeficientu teplotní délkové roztažnosti měděné trubky. Teplotní rozdíl od 20 ° C je x, prodloužení tyče je měřená veličina Y. Příklad 3 (Datový soubor palivo.txt) Pomocí (obecné) regresní paraboly spočítejte MNČ-odhady vektoru parametrů fí, aproximace Y, reziduálni součty čtverců Se a s2 v LRM (Y,X [í,cr21) pro data X 40 50 60 70 80 90 100 Y 6.1 5.8 6.0 6.5 6.8 8.1 10.0 Jedná se o měření závislosti spotřeby paliva, Y v 1/100 km motorového vozidla na rychlosti, x v km/h, při zařazeném stejném rychlostním stupni. M5120 - cvičení - Lineární regresní model 13/24 LRM - příklady (3) Příklad 4 (carb_dio. tab) Zkoumejte závislost koncentrace CO2 v atmosféře v letech 1764-1995 pomocí regresní funkce tvaru vhodného polynomu a exponeneciální funkce. Predikujte hodnoty pro 21. století. Příklad 5 (carbon_e. tab) Zkoumejte závislost uhlíkových emisí v letech 1950-1995 pomocí regresní funkce tvaru vhodného polynomu. Predikujte hodnoty pro 21. století. Příklad 6 (globtemp.txt) Zkoumejte závislost průměrné teploty v letech 1866-1996 pomocí regresní funkce tvaru vhodného polynomu. Predikujte hodnoty pro 21. století. Příklad 7 (oil_prod.txt) Zkoumejte závislost objemu vytěžené ropy v letech 1880-1988 pomocí regresní funkce tvaru vhodného polynomu a exponeneciální funkce. Predikujte hodnoty pro 21. století. M5120 - cvičení - Lineární regresní model 14/24 LRM - příklady (4) U vícerozměrných dat dále spočítejte výběrové korelační koeficienty a výběrové parciální korelační koeficienty. Obdržené hodnoty intepretujte. Příklad 8 (population.txt) Zkoumejte závislost velikosti populace na Zemi na čase pomocí regresní funkce tvaru vhodného polynomu a exponeneciální funkce. Predikujte hodnoty pro 21. století. Příklad 9 (beh.txt) Příklad na vícerozměrnou regresi (maratónské běžkyně, 1977). Pomocí regresních přímek zkoumejte závislosti mezi třemi veličinami: rychlost běhu, koroková frekvence, délka kroku. Příklad 10 (deti.txt) Příklad na vícerozměrnou regresi. Pomocí regresních přímek zkoumejte závislosti mezi třemi veličinami: hmotností dítěte, věkem a počtem bodů z diktátu. M5120 - cvičení - Lineární regresní model 15/24 LRM - příklady (5) Příklad 11 (domacnostil957.txt) Příklad na vícerozměrnou regresi (ČSR, 1957). Pomocí regresních přímek zkoumejte závislosti mezi třemi veličinami: počtem členů domácnosti, příjmy a výdaji. Příklad 12 (enrollment.txt) Příklad na vícerozměrnou regresi (VŠ v USA). Pomocí regresních přímek zkoumejte závislosti mezi veličinami: počet přihlášek na vysokou školu (ROLL), míra nezaměstnanosti (UNEM), počet absolventů střední školy (HGRAD) a průměrný příjem (INC). M5120 - cvičení - Lineární regresní model 16/24 Rozdělení pravděpodobnosti v LRM ► máme bodové odhady /S, chceme intervalové odhady a testování hypotéz ► teorie klasického LRM předpokládá e ~ N(0,tT2Jn) ► Potom máme: Y ~ N (Xj8, cr2In) ► MNČ-odhad /S vektoru parametrů je nestranný, /S ~ N (j8, cr^X'X)-1) ► s2 je nestranným odhadem rozptylu náhodných chyb, E(s2) = cr2 ► náhodné veličiny /Sas2 jsou stochasticky nezávislé ► T = c'ft-c'ft0 ^ R/ s ^'(X'X)"^ /T - j6*°Y W"1 f/T - j6*° ► F = ^-^-^-'-~ F(m, n - Jk) /S* je subvektor o m složkách W je tomuto subvektoru odpovídající blok m x m matice (X;X)_1. horní index 0 značí zvolený číselný vektor, např. při testování významnosti dosazujeme j6° = (0, • • • ,0*)', resp. j6*° = (0, • • • ,0m)' Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 17 / 24 Testy významnosti v LRM Test významnosti koeficientu j6z-, tj. test Hq: = 0 proti H\\ 7^ 0: ► v T volíme c = e\ (z-tý jednotkový vektor), j6° = 0 ► T,- = &^~t(n-fc), ^={(X'X)^},. Test významnosti modelu, Ho: (jSi,... ,jSp) = 0 proti H\: 3i > 0: j6, ^ 0: ► v F volíme j8* = (ft,...,j3p), m = k-l, j6*° = (0,...,0P) ► F= ^.|l~F(p,n-fc) =F(fc-l,n-fc) Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 18/24 Obecné testy parametrů v LRM Test lineární kombinace koeficientů Ho: c1 fS = dfS° proti H\\ c'j6 7^ c'j8° ► c G R volíme podle požadované lineární kombinace ► j8 volíme tak, aby c'jS G R byla testovaná hodnota ► T =-/ p ~ t(n - k) s y/c'[X'X)-^c ► Príklad (parabola) • B0 + B\ = 1? => volíme c = (1,1,0)',d0° = 1 ► Príklad (přímka) • 2B0 - 3B\ = 10? => volíme c = (2,3)',djf3° = 10 Vektorový test koeficientu H0: /T = j8*° proti Hx: /T / j8*°: ► testujeme subvektor j6* o m složkách (m < /c) (jp_0*oV w-i (jp_£*o\ ► F = ^-^-^-^ - F(m, n - fc) ► W je testovanému subvektoru odpovídající blok m x m matice (X'X)-1 ► Príklad • (B0,B2)' = (0,0)'? ► Príklad • (ft,,Pi)' = (0,1)'? Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 19 / 24 Korelační koeficient Korelační koeficient (Pearsonův) Rxy- _ CpCY) ► Rxy e [-1;1] ► je mírou lineární závislosti náhodných veličin X,Y (s kladnými rozptyly) ► Kovariance: C(X, Y) = E [X - E(X)] [Y - E(Y)] Pro normálně rozdělené náhodné veličiny lze interpretovat pomocí LRM s regresní funkcí Y = Čo + jSiX: ► Rxy = 1 => pozitivní lineární závislost, jSi > 0 je významný ► Rxy = —1 => negativní lineární závislost, jSi < 0 je významný ► Rxy = 0 =>- lineární nezávislost, jSi není významný Obecně platí implikace: X, Y stochasticky nezávislé ^> Rxy = 0 (nezávislost nekorelovanost) Pro normálně rozdělení náhodně veličiny platí ekvivalence (nezávislost = nekorelovanost) Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 20/24 Výběrový korelační koeficient Výběrový korelační koeficient rxy: ► je empirickou analogií korelačního koeficientu Rxy ► pro náhodný výběr ((Xi, Y\)r,..., (Xn, Yn)f) ► rXY = ^ SxSy ► Výběrový rozptyl: S2X = ^ E*=i(X ~W = lh (L/Li xí ~ "x ► Výběrová kovariance: Sxy = H=i (X - X) (Yť - Y) = ^ (E?=i X Y - nXY) ► Testy pro normálně rozdělené náhodné veličiny: ► Test nezávislosti Hq: Rxy = 0 proti H\\ Rxy 7^ 0: T = fxY Vň^2 - t(n - 2) V1_fxy ► Test Hq: Rxy = Ro proti H\\ Rxy 7^ Ro (tzv. Z-transformace): 2 1 — rxy 2 1 — Ro 2(n — 1) V-v-' z Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model Výběrový koeficient mnohonásobné korelace Výběrový koeficient mnohonásobné korelace ry.x- ► rY-x = RYXRxXRXy ► je empirickou analogií koeficientu mnohonásobné korelace Ry-x ► rY.x « Ry.x = R(Y,Y) ► Ry-x e [0; 1] ► Koeficient determinace: R2 = ľy.x ► pro náhodný výběr ((Y^Xi)7,..., (Yn,Xn)f) z (k + l)-rozměrného rozdělení ► Test pro normálně rozdělené náhodné veličiny Hn: Ry-x = 0 Pr°ti Ry-x 7^ 0» tj. že Y nezávisí na komplexu náhodných veličin X: F=-=--~ F(fc,n-fc-l) fc 1 — r2 M5120 - cvičení - Lineární regresní model Výběrový koeficient parciální korelace Výběrový koeficient parciální korelace Ty,z-x'-_ ^yz - RyxRxxRxz ► ry,z-x — / / = - RYXRX\RXY) (l - RzxRXxRxz ► je empirickou analogií parciálního korelačního koeficientu Ry7z-x ► rY/Z.x « RY/Z.X = R(Y - Y,Z - Ž) ► Ry,z-x e [-1;1] ► pro náhodný výběr ((Y^X^Zi)7,..., (yn,Xn,Zn);) z (fc + 2)-rozměrného rozdělení ► Test pro normálně rozdělené náhodné veličiny H$: Ry7z-x = 0 Pr°ti Hl: Ry,z-x 7^ 0» tj. že Y,Z jsou nezávislé náhodné veličiny po odečtení (lineárního) vlivu komplexu náhodných veličin X: T = , ľY'Z'x ~ t(n - k - 2) 2 ry,z-x Ondřej Pokora, PřF MU (2015) M5120 - cvičení - Lineární regresní model 23/24 Výběrové korelační koeficienty v R ► předpokládáme, že v proměnných X, Y, Z máme realizace náhodného výběru ► cor (X, Y) ^> výběrový korelační koeficient ľxy ► Výběrový koeficient mnohonásobné korelace ľz-{x,Y)' ► mZ <- lm (1 + X + Y) LRM na proměnných X,Y ► Zh <- mZ$fitted. values => Ž = j60 + pxX + jSyY ► cor (Z, Zh) rz.(x,Y)' = R(Z,Ž) ► Výběrový koeficient parciální korelace ľyz-x ► mY <- lm (1 + X) LRM na proměnné X ► mZ <- lm (1 + X) LRM na proměnné X ► Yr <- mY$residuals => rezidua Y - Y = Y - j60- jSxX ► Zr <- mZ$residuals rezidua Z — Z = Z — (Xq — oczX ► cor (Yr, Zr) rY/Z.x = R(Y - Y,Z - Ž) M5120 - cvičení - Lineární regresní model 24/24