MA012 Statistika II 10. Konkrétní GLM modely _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) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 1/26 Typy náhodných veličin Nominální Ordinální Intervalová Poměrová Kvalitativní Kvantitativní Dichotomická Polytomická Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 2/26 Spojitá data s normálně rozdělenými chybami Předpokládejme, že Yz- ~ N(^z-,cr ), z = 1,... ,n. Přirozeným parametrem je střední hodnota, 0Z- = ^. Pro modelování střední hodnoty EYZ = ^ používáme identickou linkovací funkci g = id, g(m) = m = 9i = tu = ýp. Takový GLM vede na klasický lineární regresní model. Gaussova-Markovova věta zaručuje, že odhady parametrů j6 metodou nejmenších čtverců jsou BLUE (nejlepší lineární nestranný odhad, best unbiased linear estimator), tedy pro LRM s parametry /S je dosažen nejmenší reziduálni součet čtverců. Normalita dat zaručuje, že odhady j6 metodou nejmenších čtverců jsou stejné, jako odhady j6ML metodou maximální věrohodnosti. MA012 Statistika II - 10. GLM - Konkrétní GLM modely 3/26 m Alternativní a binomická data Předpokládejme, že li; ~ A(pi), i = 1,...,N, nabývá pouze dvou hodnot 0 a 1, Vi P(Ui = u)={ 1-p, 0 u = 1 M = 0 = jinak Pľ(l-Pi)1_w " = 0,1 o jinak. Předpokládejme, že náhodná veličina Uj závisí na k kovariátech Xi\,... Data můžeme mít zadána různým způsobem: jednotlivá pozorování Uf hodnoty kovariát pozorované binární veličiny Xj\, . . ., Xjfc Ui skupinově - pomocí absolutních četností úspěchů Yy v rtj pokusech -,a- J Q)pí i1 - Pi)nry y = 0'1..... P(Yj = y) = n j 0 jinak kde j = 1,... ,n, N = n\ + • • • + nn, a data můžeme zapsat formou tabulky hodnota kovariát počet úspěchů počet pokusů Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 4/26 Alternativní a binomická data skupinově - pomocí relativních četností úspěchů Z; = v m pokusech [ 0 jinak kde j = 1,... ,n, N = ni + • • • + nn, a data lze zapsat do tabulky kovariáty relativní úspěšnost počet pokusů Z j Yj/tij Úkolem statistické analýzy je nalézt vztah mezi Zj (resp. Y f) a kovariátami V GLM modelujeme pravděpodobnosti pj pomocí linkovací funkce g(Pi) = Vi = xiP-Nejjednodušším modelem je lineární model s g = id, Pi = *íj6. Avšak tento model má řadu nevýhod, především je třeba zajistit, aby x-/S nabývala hodnot mezi 0 a 1, tedy je třeba přidat nějaké dodatečné podmínky. Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 5/26 Probitový model Probitový model používá tzv. probitovou linkovací funkci g, což je kvantilová funkce O-1 standardizovaného normálního rozdělení. GLM je tedy ve tvaru g(pi) = O^ipi) = uPi = rji = x[fi. To znamená, že pravděpodobnost úspěchu je modelována vztahem kde Oje distribuční funkce standardizovaného normálního rozdělení. V případě rji = /3o + /3i x\ potom dostáváme pz- = (/3o + $\ x\) = F(xj), kde F je distribuční funkce rozdělení N(—/3o//3i, jS^2). f(s)~ N^o2) jc(x) ~ Ním-.o2) Logistický (logitový) model - logistická regrese Logistický model používá tzv. logitovou linkovací funkci g, což je kvantilová funkce logistického rozdělení pravděpodobnosti. GLM je ve tvaru g(Pi) = In j^— = rji = x[p. 1 yi Pravděpodobnost úspěchu, resp. neúspěchu, je modelována vztahem vi = g \ví) = y + exp(-x/j6)/ 1 ~Pi l + exp(x7j6) (i (i Obr. Porovnání funkcí g_1 v probitovém (plně) a logitovém (čárkovaně) GLM Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 7/26 C Log Log (komplementární LogLog) model CLogLog model používá jako linkovací funkci g kvantilovou funkci logaritmického Weibullova (extreme-minimal-value) rozdělení pravděpodobnosti. GLM je ve tvaru g(Pi) = ln [- ln(! - Vi)\ = Vi = xiP- Pravděpodobnost úspěchu je tedy modelována vztahem Vi = S 1(Vi) = 1 - exp [- exp(^jS) 7t(x) Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 8/26 LogLog model LogLog model používá jako I in kovací funkci g kvantilovou funkci zobecného Gumbelova (extreme-maximal-value) rozdělení pravděpodobnosti. GLM je ve tvaru g(Pi) = - ln (- lnp) = jfc = xiP-Pravděpodobnost úspěchu je tedy modelována vztahem Vi = S 1(Vi) = exP [- exp(-^jS) Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 9/26 Logistická regrese a šance Definice 1 (Šance) Podíl pravděpodobnosti úspěchu a pravděpodobnosti neúspěchu v binomickém (příp. alternativním) rozdělení pravděpodobnosti se nazývá šance (odds), odds = P(Y = 1) p P(Y = 0) 1-p Logaritmus ln odds je v finančnictví někdy nazýván skóre. V logistické regresi se používá logitová linkovací funkce g(p) = lny^ Pravděpodobnost úspěchu, resp. neúspěchu, je pak rovna P g W l + exp^)' P = l + exp(*'j6)' a šance vychází jako exponenciála lineárního prediktoru odds = , 6XP(X^L = exp(^) = e", ln «fcfe = 7 l + exp(x'£) 1 ťV p; ' = r/ = x' p. MA012 Statistika II - 10. GLM - Konkrétní GLM modely Logistická regrese a podíl šancí Pokud v logistické regresí modelujeme binární veličinu Y £ {0, 1} pomocí taktéž binární kovariáty x £ {0, 1}, měříme asociovanost těchto dvou veličin pomocí statistiky QR: Definice 2 (Podíl šancí) Podíl šancí (odds ratio, podíl rizik) je podílem podmíněných šancí, or = odds(x 1) odds(x 0) P(Y= =1 x= =1) P(Y= =0 x= =1) P(Y= =1 x= =0) P(Y= =0 X- =0) P(Y = 1 x = l)P(Y = 0 x = 0) P(Y = 0 x = 1) p(Y = 1 x = 0) OR « 1: nekorelovanost OR > 1: korelovanost, souhlasná asociovanost OR < 1: korelovanost, nesouhlasná asociovanost Dosazením šancí dostáváme výpočetní tvar odds(x = 1) exp(V/S) x=l odds(x = 0) exp(x7j6) | x=q Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 11/26 Logistická regrese a podíl šancí V případě modelu logistické regrese veličiny Y pomocí binární proměnné x G {0, 1} s lineárním prediktorem tvaru r] = jSq + j^i x, tzn. v GLM Y ~ Bi{n, p), ln v v dostáváme or = exp(j60 +j6i*) x=l exp(j60 +j6i*) |x=0 exp(j80 + j8i) e*p(/30) = eči tedy lineární koeficient /3i binární kovariáty x v modelu logistické regrese je rovný logaritmu podílu šancí, /3i = ln QR. MA012 Statistika II - 10. GLM - Konkrétní GLM modely Poissonovská data Předpokládejme, že náhodný výběr rozsahu n je z Poissonova rozdělení, tj přičemž EYZ = DYZ = Az-. Yj = počet výskytů sledovaného jevu v určitém časovém intervalu, na dané ploše, v daném objemu, apod. Podmínky: a) jev může nastat v kterémkoliv časovém okamžiku, )) počet výskytů jevu během časového intervalu závisí jen na jeho délce a ne na jeho počátku ani na tom, kolikrát jev nastoupil před jeho počátkem, c) pravděpodobnost, že jev nastoupí více než jednou v intervalu délky t, konverguje k nule rychleji než t, d) A je střední hodnota počtu výskytů jevu za časovou jednotku. \i > 0; y = 0,1,2,... jinak MA012 Statistika II - 10. GLM - Konkrétní GLM modely Modely pro Poissonovská data Logaritmicko-lineární model používa logaritmickou linkovací funkci g. GLM je ve tvaru g(\i) =ln\i = ?ii = x/ip. Střední počet udalostí je tedy modelován exponenciálně, Odmocninový model používá lin kovací funkci g ve tvaru odmocniny, GLM je ve tvaru Střední počet událostí je tedy modelován kvadraticky K = g-\ni) = {x'iP)2- MA012 Statistika II - 10. GLM - Konkrétní GLM modely Další modely Pro data s exponenciálním a gama rozdělením pravděpodobnosti, tj. pro spojitá data s nekonstantní chybou a konstantním poměrem střední hodnoty a směrodatné odchylky, se obvykle používá inverzní linkovací funkce g. GLM je ve tvaru 1 g(m) = — = Vi = Ač-ri Střední hodnota je tedy modelována vztahem Pro data s inverzním Gaussovým rozdělením pravděpodobnosti se obvykle používá inverzní kvadratická linkovací funkce g. GLM je ve tvaru 1 Střední hodnota je tedy modelována vztahem m = g~1(Vi) = Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 15/26 Zobecněné lineární modely v R Obecná funkce pro řešení GLM v r je glm. model <-glm (formula, family, data) family family (link = ...) gaussian identity, log, inverse binomial logit, probit, cloglog, log, cauchit poisson log, sqrt, identity Gamma inverse, log, identity inverse.gaussian l/mu~2, inverse, log, identity v <- summary (model) S výsledky se pracuje analogicky jako s výsledky funkce lm pro LRM. MA012 Statistika II - 10. GLM - Konkrétní GLM modely Příklad Příklad 1 V souboru beetle.csv jsou uvedeny údaje o úmrtnosti Potemníka skladištního (Tribolium confusum) v reakci na sirouhlík CS2. Datový soubor obsahuje tyto proměnné do se množství sirou hlíku (mg/l) population počet kusů ve zkoumaném vzorku kil led počet mrtvých kusů ve zkoumaném vzorku Modelujte závislost úmrtnosti na množstvíCS^- Řešení. Pro modelování závislosti použijeme logistický model, probitový model a model s komplementární log-log linkovací funkcí. MA012 Statistika II - 10. GLM - Konkrétní GLM modely 17/26 Link Function g(jt)=log(jt/(1-Jt)) Normal Probability Plot y=*(-35.127+19.838x) link function g(jt)=log(-log(1-Jt)) 1.7 1.75 1.8 1.85 y=1 -exp(-exp(-40.647+22.656x)) -3-2-1 0 1 2 3 4 Standard Normal Deviate Normal Probability Plot K 0 -3-2-1 0 1 2 3 4 Standard Normal Deviate Normal Probability Plot K 0 -3-2-1 0 1 2 3 4 Standard Normal Deviate ik: Modely pro úmrtnost Potemníka skladištního. Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 18/26 Příklad Příklad 2 V souboru aids. csv jsou uvedeny údaje o počtech nových případů AIDS ve Velké Británii za období prosinec 1982 až listopad 1985. Datový soubor obsahuje tyto proměnné month měsíc year rok nurnber počet nových případů AIDS Modelujte závislost počtu nových případů AIDS na čase. Řešení. Pro modelování závislosti použijeme lineární model, log-lineární model a odmocninový model. MA012 Statistika II - 10. GLM - Konkrétní GLM modely 19/26 link function g(n)=n Normal Probability Plot • • • 10 15 20 25 30 35 y=-2.2+0.457x link function g(n)=log(n) 5 10 15 20 25 30 35 y=exp(0.0397+0.0796x) link function g(n)=Vn -1-1-1-1-w~r -3-2-10123 Standard Normal Deviate Normal Probability Plot -3-2-10123 Standard Normal Deviate Normal Probability Plot -1-1-1- —i-1-1-1- y ++^ y y y y y y=(0.553+0.0944xr -3-2-1 0 1 2 Standard Normal Deviate ik: Modely pro výskyt nových onemocnění AIDS ve Velké Británii Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétni GLM modely 20/26 Neškálová deviance Předpokládáme, že náhodný výběr Y = (Y\,...,Yn) se řídí GLM modelem, t i=l i=l Předpokládejme dále, že hustota je ve škálové formě, tj. %{ 0, CVj se známými váhami coj > 0 a jedním neznámým rušivým parametrem (scale factor)

0. Škálová deviance je rovna D = 2 £(L;Y)-£(Č;Y) z=l max #i) - l(6i,max) + 7(^i)] = D* d* a D* nazýváme neškálovou deviancí (unsealed deviance). Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 21/26 Odhady rušivého parametru Protože platí D* as. 2/ n _^ t-t-» ED* D = — ~ r n — /c ED =-« n lze rušivý parametr ^ odhadnout pomocí statistiky D* <£d* = -r- n —k Pro tzv. zobecněnou Pearsonovu statistiku x2 platí n /v. Odtud lze získat jiný odhad rušivého parametru, x2 Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely Overdispersion, underdispersion Připomeňme, že v normálním rozdělení je rušivým parametrem rozptyl, (p = a a v gama rozdělení (p = l/k. V prostředí R je pro řešení problémů s neadekvátně odhadnutým rozptylem k dispozici modifikovaná volba pro třídu exponenciálního rozdělení. V případě binomického rozdělení máme možnost volby family = quasibinomial a pro Poissonovo rozdělení family = quasipoisson. Nejedná se přiotm o nové rozdělení exponenciálního typu, ale o změnu ve výpočtu druhého momentu, pro jehož odhad se použije jednoduchý momentový odhad disperzního parametru (p. Výsledná korekce rozptylu je pak důležitá při testování hypotéz, neboť zohledňuje vyšší či nižší variabilitu v datech a zabraňuje tak nadbytku či nedostatku falešně pozitivních výsledků testů hypotéz o parametrech modelu. MA012 Statistika II - 10. GLM - Konkrétní GLM modely Příklad Příklad 3 V souboru bees.csv jsou uvedeny údaje o aktivitě včel v závislosti na čase. Jednou z důležitých charakteristik při zkoumání včelí aktivity je počet včel, které opustí úl kvůli práci ve vnějším prostředí. Studie se zabývala měřením této veličiny během několika slunečných dní v závislosti na čase během dne. Datový soubor obsahuje tyto proměnné number počet včel, které opustily úl time čas, kdy byl tento údaj zaznamenán Modelujte závislost počtu včel, které opustí úl, na čase během dne. Řešení. Pro modelování závislosti použijeme poissonovský model s kanonickou I in kovací funkcí. Do modelu vstupuje jediná vysvětlující proměnná time a přidáme také její druhou mocninu. Hodnota reziduálni deviance (4 879,3) je nepoměrně vyšší než počet stupňů volnosti (501). Je zřejmé, že došlo k „overdispersion" a v jazyce R je třeba volit f amily=quasipoisson. Použití této volby neovlivňuje odhady koeficientů, ale mění jejich odhady variability, což se projeví např. v intervalu spolehlivosti. Ondřej Pokora, PřF MU (2015) MA012 Statistika II - 10. GLM - Konkrétní GLM modely 24 / 26 Bees activity 8 10 12 14 16 time ik: Odhad regresní funkce bez vyrovnání se s problematikou velkého rozptylu. MA012 Statistika II - 10. GLM - Konkrétní GLM modely Bees activity 8 10 12 14 16 time ik: Odhad regresní funkce s vyrovnáním se s problematikou velkého rozptylu. MA012 Statistika II - 10. GLM - Konkrétní GLM modely