M5VM05 Statistické modelování 10. Konkrétní GLM modely - I.
Jan Koláček (kolacek@math.muni.cz)
Ustav matematiky a statistiky, Přírodovědecká fakulta, Masarykova univerzita, Brno
Jan Koláček (PřF MU) M5VM05 Statistické modelování 1/33
Motivace
Na minulé přednášce jsme si uvedli obecnou definici zobecněného lineárního modelu a obecné konstrukce testů hypotéz o parametrech těchto modelů. Na této přednášce se již budeme zabývat zobecněnými lineárními modely pro konkrétní případy podle toho, jaké rozdělení má závisle proměnná Y.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
2/33
Motivace
Typy veličin:
Nominální Ordinální
Kvalitativní
Intervalová Poměrová
Kvantitativní
Diskrétní 2
Spojitá
Kategoriální
D ichot omická     Poly t omická
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
Modely pro alternativní a binomická data
Předpokládejme, že Uj ~ A(tzí) (i = 1,... ,N) nabývá pouze dvou hodnot 0 a 1,
nf(l - tzí)1-11   u = 0,1
7T; U = 1
P(Ui = «)=<(    1 - 7Tŕ    U = 0
0 jinak
0
jinak.
Předpokládejme, že náhodná veličina Uj závisí na A: veličinách xn,... ,x^, tzv. kovariáty.
Data můžeme mít zadána různým způsobem: • jednotlivá pozorování U,\
hodnoty kovariát	pozorované binární veličiny
%il i ■ ■ ■ i %ik	Ui
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
Modely pro alternativní a binomická data
o skupinově, kdy známe absolutní četnosti úspěchů Yj a celkový počet pokusů rij, tedy máme k dispozici binomická data Yj ~ Bi(rij,TZj)
PiY1=y)={ W V-y = M"/ K 1    i"     \ 0 jinak
kde
j = 1,..., n;      N = n1-\-----h n„
a data můžeme zapsat formou tabulky
hodnota kovariát	počet úspěchů	počet pokusů
Xji, . . ., Xjfc		rij
Jan KoláCek (PřF MU)
M5VM05 Statistické modeloval
5/33
Modely pro alternativní a binomická data
• skupinově, kdy známe relativní četnost úspěchů Zj =     a celkový počet pokusů tij
q^il-n^   y=0,i.....1
0 jinak
kde
j = 1,..., n;      N = n1-\-----h n„
Data lze zapsat do tabulky
kovariáty	relativní úspěšnost	počet pokusů
Xji,. . .,Xjfc	zř = £ 1 tij	rij
P(Zj = y) =
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
6/33
Cíl
Hlavním úkolem statistické analýzy je nalézt vztah mezi Z,-, (tj. i Y,) a xn,.. . ,x,7t, tj. funkci
7T; = 7r(Xi) = 7T(xfl,...,Xr;t).
Protože chceme použít GLM modely, modelujeme pravděpodobnosti 7i(- pomocí linkovacích funkcí
Nejjednodušším modelem je lineární model
Tli = Xjj8.
Avšak tento model má řadu nevýhod, především je třeba zajistit, aby xj/3 nabývala hodnot mezi 0 a 1, tedy je třeba přidat nějaké dodatečné podmínky. Proto, abychom tuto podmínku dodrželi, využijeme nějakou distribuční funkci
f{s)ds     f(s) > 0       /   f(s)ds = 1
-co J —co
s odpovídající hustotou f(s), která se v tomto případě nazývá toleranční funkce (toleranční distribuce).
Jan Koláček (PřF MU) M5VM05 Statistické modelování 7/33
Modely dávka - odpověď
Typickým příkladem těchto modelu je vztah mezi dávkou toxické látky a odezvy (kladná-přežití, záporná-smrt) jedince na tuto dávku. Odezvy bývají obvykle udávány jako procenta kladné odezvy (quantal responses). Symetrické modely
Jestliže uvažujeme toleranční distribuci jako rovnoměrně spojitou na nějakém intervalu (a,b), tj
pak
pro   x e (a, b)
a tento model je lineárním modelem
7T0(x)
X —
a
Po + fax
tj- ft,
a
1
b — a
b — a
b — a
> 0
s identickou linkovací funkcí
^o(tt) = n.
Jan KoláCek (PřF MU)
M5VM05 Statistické modelovaní
8/33
Symetrické modely
Obrázek : Rovnoměrné rozdělení na (a,b).
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
9/33
Probitový model
Další možností je vzít normální hustotu jako toleranční funkci. V tomto případě mluvíme o tzv. probitovém modelu:
nliX) = FlW = £/l(s)ds = £ ^e-H^ds = * (^),
kde Oje distribuční funkce standardizovaného normálního rozdělení. Pak tzv. probitovou linkovací funkcí je kvantilová funkce normálního rozdělení
gl(n) = <S> 1(n)
x— cr
Čo + jM     tj. 0O
h = l > o.
Hodnota mediánu x = \i se nazývá mediánová smrtící dávka (median lethal dose - LD50) a odpovídá dávce, při které polovina jedinců má kladnou a polovina zápornou odezvu.
Jan Koláček (PřF MU)
M5VM05 Statistické modelováni
Probitový model
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
11 / 33
Logistický model
Jiným velmi podobným modelem je tzv. logistický model, kde toleranční funkce je hustota logistického rozdělení
f (s) = I_ÍÍPÍ5lÍ)_ = i ^pQ^it) nK>      -[l+exp(^)]2      - [l+exp(-^)]2'
takže
7T2(x) = F2(x)
i    eMs-ir)    ds =   exP(V)   = i
s tzv. logit linkovací funkcí
g2(7i)=log(I^) = ^ = /30 + /31x      tj.      /$0 = -£   /31 = i>0.
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
Logistický model
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
13 / 33
CLogLog model
Asymetrické (extremální) modely
Pokud za toleranční funkci zvolíme Log-Weibullovo rozdělení
(extreme-minimal-value distribution) ve tvaru
/3(s) = i exp (^/J exp
pak
7T3(x) =
s tzv. komplementární log-log linkovací funkc
exP \a
iexpí^/lexp
■ exp ^
ds = 1 — exp
exp
g3(7r) = log[-log(l-7i)] = Vi = ßo + /3iX   tj.   ß0 = -%   /31 = i>0.
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
CLogLog model
Jan Koláček (PřF MU)
M5VM05 Statistické modelová
LogLog model
Pokud jako toleranční funkci zvolíme zobecněné Gumbelovo rozdělení
(extreme-maximal-value distribution) ve tvaru
/4(s) = i exp exP
dostaneme
s tzv. log-log linkovací funkcí
g3{n) = -log[-log(7r)] =
exp
cr
i exp
cr
exp
s—]l cr
ds = exp
x— cr
Čo + jM tj.
exp
-l   č1 = I>0.
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
LogLog model
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
Logistická regrese
Nejčastěji se používá logit linkovací funkce
g2(7l) = log(T^).
Zajímá nás vztah pravděpodobností úspěchu či neúspěchu k hodnotám regresorů (kovariát) x= {x\,... ,X;t)T, tj.
P(Y-l\xi       xt)-n(x)-    exP^(x)>    -_l_
111*k)-nW- 1+eXp{7/(x)} ~ l+exp{-?7(x)}
a
P<Y=°i-.....*> =1 - *<*> = Tr4mi= tw-,^)}
Předpokládejme, že lineární prediktor je roven
?/(x) = /30 + /3Tx.
Jan Koláček (PřF MU) M5VM05 Statistické modelování
Logistická regrese
Všimněme se nejprve, že podíl
odds(l) _ P(Y = l\xlr...,xk) _ 7i(x)
odds(O)     P(Y = 0\x1,...,xk) l-7i(x)
exp(/30 + /3Tx)
má bezprostřední interpretaci. Porovnává pravděpodobnost jedničky (tj. výskyt sledovaného jevu při daných hodnotách kovariát) a nuly (nevýskyt sledovaného jevu při daných hodnotách kovariát). Anglickému označení odds odpovídá české označení šance. Pro k = 1 jsou šance
odds(O) = exp(/3o), odds(l) = exp^o + jSi)-Poměr šancí (anglicky odds ratio) pro binární x je pak
^„    odds(í) ,n .
takže parametr f$i je roven logaritmu poměru šancí.
Jan KoláCek (PřF MU)
M5VM05 Statistické modelováni
Příklad
Příklad 1
V souboru „beetle.RData" jsou uvedeny údaje o úmrtnosti Potemníka skladištního (Tribolium confusum) v reakci na sírouhlík CS2. Datový soubor obsahuje tyto proměnné
dose   množství sírouhlíku (mg/l) population   počet kusů ve zkoumaném vzorku
killed   počet mrtvých kusů ve zkoumaném vzorku
Modelujte závislost úmrtnosti na množství CSi-
Řešení. Pro modelování závislosti použijeme logistický model, probitový model a model s komplementární log-log linkovací funkcí.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
20 / 33
Příklad
Obrázek : Modely pro úmrtnost Potemníka skladištního.
Jan Koláček (PřF MU) M5VM05 Statistické modelování
21 / 33
Modely pro poissonovská data
Předpokládejme, že náhodný výběr rozsahu n je z Poissonova rozdělení, tj
Y, ... počet výskytu sledovaného jevu v určitém časovém intervalu (na ploše velikosti t apod.).
Jestliže jsou splněny následující podmínky
a) jev může nastat v kterémkoliv časovém okamžiku,
b) 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
pak uvedená náhodná veličina má rozdělení Po(A).
Y,-~Po(A,-),   P{Yi = y)
0
A,- > 0; y = 0,1,2,... jinak
přičemž
EYi = DYj = A,-.
Jan KoláCek (PřF MU)
M5VM05 Statistické modeloval
22 / 33
Modely pro poissonovská data
Náhodnou veličinou, která má Poissonovo rozdělení, je tedy např.
• počet vadných výrobku ve velké sérii, jestliže pravděpodobnost vyrobení vadného výrobku je velmi malá
• počet těžkých dopravních úrazů za den v určitém městě
• počet zákazníků v prodejně během nějakého časového intervalu
• počet částic v jednotce plochy nebo objemu, např. počet částic v zorném poli mikroskopu
o počet telefonních volání v časovém intervalu t
• počet létavic pozorovaných během intervalu délky t
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
23 / 33
Modely pro poissonovská data
Předpokládejme, že Y,- závisí na k veličinách xn,... ,x,7t, a úkolem je najít vztah mezi nimi, tj. hledáme funkci
A,- = A(x,-) = A(xilr...,xik).
GLM modely =>■ modelujeme pravděpodobnosti A,- pomocí linkovacích funkcí
Nejjednodušším je lineární model
A,- = xTjg.
Tento model má řadu nevýhod, především je třeba zajistit, aby nabývala pouze kladných hodnot. Nejčastěji se volí tyto dvě možnosti: log-lineárních model :
EYi = fa = A{ = exp(xJ/3) (#) = rji = gx (A,) = log(A,) =
odmocninový model :
EY; =      = A; = (x]fi)2 g2{}li) = f]i = g2{Ai) = ^Ä" =
Jan KoláCek (PřF MU)
M5VM05 Statistické modelování
24
Modelování binomických dat pomocí poissonovského modelu
Pomocí Poissonova rozdělení Po(A) lze dobře aproximovat binomické rozdělení Bi(n, 7i) za podmínek
n —>■ oo   & 71—^0   &   nn —>■ A < oo,
obvykle se doporučuje n > 30 a n < 0,1.
Chceme-li tedy aproximovat binomické rozdělení Bz(n,-, 7T() pomocí Poissonova rozdělení Po(A(- = n,-7T() a přitom použijeme logaritmickou linkovací funkci, platí
A,- = n,-7T,- = exp(xf/S)     log(A,) =   log(nO  + log(7rr) = xf/S.
tzv. „offset"
Jan KoláCek (PřF MU)
M5VM05 Statistické modeloval
25 / 33
Příklad
Příklad 2
V souboru „aids.RData" 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 number   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.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
26 / 33
Příklad
Obrázek : Modely pro výskyt nových onemocnění AIDS ve Velké Británii.
Jan Koláček (PřF MU) M5VM05 Statistické modelování
27 / 33
Úlohy k procvičení
Příklad 3.1
V souboru „heart .RData" jsou uvedena data o přítomnosti infarktu myokardu v závislosti na věku pacienta. Datový soubor obsahuje tyto proměnné:
age    věk pacienta (roky)
chd   indikátor infarktu (1 - nastal, 0 - nenastal)
Pro modelování závislosti použijte logistický model, probitový model a model s komplementární log-log linkovací funkcí. Výsledky vykreslete do obrázku.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
28 / 33
Úlohy k procvičení
Příklad 3.2
V souboru „nemocnice.RData" jsou uvedeny údaje o zotavení pacientů v závislosti na závažnosti onemocnění a nemocnici, ve které se léčili. Datový soubor obsahuje tyto proměnné:
InfectionSeverity    vážnost onemocnění Treatment_Outcome     indikátor uzdravení (1 - zdravý, 0 - smrt) Hospital typ nemocnice (1, 2, 3)
Pro modelování závislosti nalezněte vhodný logistický model. Výsledky vykreslete do obrázku.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
29 / 33
Úlohy k procvičení
Příklad 3.3
V souboru „cancer.RData" jsou uvedeny údaje o počtu onemocnění rakovinou kůže u žen v závislosti na věku a oblasti v USA, ve které pacientky žily. Datový soubor obsahuje tyto proměnné:
Cases počet onemocnění
Town město (0 - Minneapolis (Minnesota), 1 - Dallas (Texas))
Age věková skupina pacientky
Population celkový počet žen dané věkové skupiny v příslušném městě
Pro modelování závislosti nalezněte vhodný logistický model. Výsledky vykreslete do obrázku. Porovnejte pravděpodobnost vzniku onemocnění u 60-ti leté pacientky žijící v Minneapolisu s pravděpodobností pro stejně starou pacientku žijící v Dallasu.
[Minneapolis: 0.00117, Dallas: 0.00276/
Jan Koláček (PřF MU) M5VM05 Statistické modelováni
Úlohy k procvičení
Příklad 3.4
V souboru „ car-income. RData" jsou uvedeny údaje o koupi nového auta během posledních 12-ti měsíců v závislosti na přijmu domácnosti a stáří původního auta. Datový soubor obsahuje tyto proměnné:
purchase   indikátor nákupu nového auta (1 - ano, 0 - ne) income       roční příjem domácnosti (v tis. dolarů) age stáří původního auta (roky)
Nejprve vykreslete závislosti proměnné purchase na ostatních. Pro modelování závislosti nalezněte vhodný logistický model. Jsou všechny proměnné statisticky významné? Znovu modelujte s použitím proměnné age jako factor. Opět sledujte statistickou významnost age. Vyzkoušejte tuto proměnnou zakomponovat do modelu jako factor s méně úrovněmi. Výsledky vykreslete do obrázku.
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
31 / 33
Úlohy k procvičení I
Příklad 3.5
V souboru „ druhy .RData" jsou k dispozici data, která se týkají dlouhodobého zemědělského experimentu. Bylo sledováno 90 pozemků (pastvin) o rozloze 25 m x 25m, lišících se v biomase, pH půdy a druhové bohatosti (počet rostlinných druhů na celém pozemku). Je dobře známo, že s rostoucí biomasou dochází k poklesu druhové bohatosti. Ale zůstává otázka, zda rychlost poklesu nesouvisí s úrovní pH v půdě. Proto byly jednotlivé pozemky klasifikovány podle hodnoty pH v půdě do tří úrovní (nízká, střední a vysoká úroveň) a do experimentu bylo vybráno vždy po 30 pozemcích pro každou úroveň. Spojitá veličina Biomass je dlouhodobým průměrem naměřených červnových hodnot biomasy. Datový soubor obsahuje tyto proměnné:
pH úroveň pH v půdě (low - nízká, mi d - střední, high - vysoká)
Biomass    množství biomasy species   počet rostlinných druhů
Jan Koláček (PřF MU) M5VM05 Statistické modelováni
Úlohy k procvičení II
Příklad 3.5
Nejprve vykreslete závislosti proměnné species na ostatních. Pro modelování závislosti nalezněte vhodný poissonovský model. Vyzkoušejte postupně logaritmickou, identickou a odmocninovou linkovací funkci. Jsou všechny proměnné statisticky významné? Pokud ne, zkuste modely zjednodušit a pomocí analýzy deviace rozhodněte, zda takové zjednodušení je možné. Získané výsledné modely vykreslete do obrázku. Pomocí všech modelů odhadněte počet rostlinných druhů na pozemku s hodnotou biomasy 9 a střední úrovní pH v půdě.
[Odhady počtu druhů pro log link: 8,895, identity link: 4,513, sqrt link: 7,414.]
Jan Koláček (PřF MU)
M5VM05 Statistické modeloval
33 / 33