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é modelování
2/33
Motivace
Typy veličin
Nominální		Ordinální	
Intervalová Poměrová
Kvalitativní
Kvantitativní
Diskrétní
Spojitá
Kategoriální
Dichotomická
Polytomická
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
3/
Modely pro alternativní a binomická data
Předpokládejme, že Uj ~ A(tZj) (z = 1,.. .,N) nabývá pouze dvou hodnot 0 a 1,
P(Lřf = u)
Tli 1
0
7Ty
U = 1 M = 0
jinak
_   í   TďOL-TTf)1-"    M = 0,1
0
jinak.
Předpokládejme, že náhodná veličina líz- závisí na k veličinách X\\,... ,Xfc, tzv. kovariáty.
Data můžeme mít zadána různým způsobem: 9 jednotlivá pozorování Uf
hodnoty kovariát	pozorované binární veličiny
X\\, . . ., Xfc	Lřf-
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
4/33
Modely pro alternativní a binomická data
skupinově, kdy známe absolutní četnosti úspěchů Yy a celkový počet pokusů íij, tedy máme k dispozici binomická data Yy ~ Bi(rij,7Zj)
P{Yj = y)
Q>/(1 - JZjp-y   y = 0,1.....n
J
0
jinak
kde
; = 1,.. .,n;      N = ni H-----hnn
a data můžeme zapsat formou tabulky
hodnota kovariát	počet úspěchů	počet pokusů
		
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
5/33
Modely pro alternativní a binomická data
• skupinově, kdy známe relativní četnost úspěchů Z, = -r^ a celkový
J ilj
pokusů rij
7 [0 jinak
kde
; = 1,.. .,n;      N = ni H-----hnn
Data lze zapsat do tabulky
kovariáty	relativní úspěšnost	počet pokusů
	] rij	tij
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
Cíl
Hlavním úkolem statistické analýzy je nalézt vztah mezi Zz-, (tj. i Yf) a X{\,... ,Xfc, tj. funkci
Tli = 7ľ(Xj) = 7ľ(Xji, . .
Protože chceme použít GLM modely, modelujeme pravděpodobnosti ttz- pomocí linkovacích funkcí
g(7li) =
Nejjednodušším modelem je lineární model
Tli = xJ'jS.
Avšak tento model má řadu nevýhod, především je třeba zajistit, aby xjfi 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
/l ľ co
f(s)ds     f(s)>0       /   /(s)ds = l -00 J CO
s odpovídající hustotou /(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 - odpoveď
Typickým příkladem těchto modelů 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
s G (a, b) jinak
pak
pro   x <G (a, b)
a tento model je lineárním modelem
7T0(x)
x — a
jSo + fr*
tj- čo
a
> o
b — a
b — a
b — a
s identickou linkovací funkcí
g0(n) = ti.
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
8/
Symetrické modely
f(s)~ Rs(a,b)
7t(x) ~ Rs(a,b)
1/(b-a)
a (a+b)/2 b
a (a+b)/2 b
Obrázek : Rovnoměrné rozdělení na (a,b)
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
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:
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)=^-\n) = x-^=^ + ^x     tj.      h = ~l   jSi = i > 0.
Hodnota mediánu x — ]i se nazývá mediánová smrtící dávka (medián 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ání
10 / 33
Obrázek : Normální rozdělení N[}L,a).
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
11 / 33
Logistický model
Jiným velmi podobným modelem je tzv. logistický model, kde toleranční funkce je hustota logistického rozdělení
[l+exp(^i)
l+exp(-^)
2/
takže
7T2(x) = F2(X) =  ľ   1 2 rfs =    6XP^J, = -1
2W       2V ;     7-oo <Ml+exp(^)]2 l+expCV1) !+exP(
s tzv. logit linkovací funkcí
^2(/r)=log(T^)
0"
0o + ft*     tj. =        ]6i = ^>0
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
12 / 33
Obrázek : Srovnání probitového a logistického (---) modelu při stejných parametrech ]i
a a.
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
13 / 33
C Log Log model
Asymetrické (extremální) modely
Pokud za toleranční funkci zvolíme Log-Weibullovo rozdělení
(extreme-minimal-value distribution) ve tvaru
/3(s)
\ exp (V1) exP
s—u
exPi V1
pak
7T3(x)
4X
00
^expf^ JexP
exp
S — fl
<J
ds = 1 — exp
exp
<J
s tzv. komplementární log-log linkovací funkcí
g3(n) =log[-log(l-7r);
_ x—]i _       O    \   o        +■      o _ V-     o _ 1
(7
po + fax  tj.   |6o = -£   jSi = ^ > O
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
C Log Log model
f(S) Jt(x)
Obrázek : Log-Weibullovo rozdělení.
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
15/
LogLog model
Pokud jako toleranční funkci zvolíme zobecněné Gumbelovo rozdělení
(extreme-maximal-value distribution) ve tvaru
/4(s) = i exp (-s-^ j exp
exp
S — fl
<J
dostaneme
TC±(x)
-00
i exp
S — fl <J
exp
S — fl
<J
ds exp
exp
X — fl
<J
s tzv. log-log linkovací funkcí
g3(n) = -log[-log(7r)
_ X fl _ a    i  n. v a„ —     V-     o   _ 1
0"
Čo + frx  tj.   & =        ]6i = ^>0
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
16/
Obrázek : Zobecněné Gumbelovo rozdělení.
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
17 / 33
Logistická regrese
Nejčastěji se používá logit linkovací funkce
ft(*)=l0g(i^ř).
Zajímá nás vztah pravděpodobností úspěchu či neúspěchu k hodnotám regresorů (kovariát) x = {x\,... ,x^)T, tj.
p(y = i \x\,..., Xfc) = tt(x)
exp{7/(x)}    _ 1 1 + exp{7/(x)}     1 + exp{—f](x)}
P(Y    0\xi,...      = 1 — tt(x) —
1 exp{—//(x)}
+ exp{7/(x)}     1 + exp{-7/(x)}
Předpokládejme, že lineární prediktor je roven
jy (x) = j80 + j6Tx
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
18 / 33
Logistická regrese
Všimněme se nejprve, že podíl
jj  í \        P(Y = 1 • -rXk) 7ľ(X) /n /,T \
0íMs(xiXfc) = Pir^ox!.....Xk) = T^W) = exp(^°+ * x)
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(0) = čxp(/3o), odds(l) = exp(Po + j8i). Poměr šancí (anglicky oc/c/s rař/o) pro binární x je pak
takže parametr /5i je roven logaritmu poměru šancí.
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
19 / 33
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é
do s e   množství sírouhlí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í CS2-
Ř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é modelování
20 / 33
Příklad
Link Function g(jt)=log(jt/(1-jt))
Normal Probability Plot
y=*(-35.127+19.838x) link function g(jt)=log(-log(1-:t))
-3-2-1 0 1 2 3 4
Standard Normal Deviate
Normal Probability Plot
++-M-
4- +
—i_i_
Normal Probability Plot
y=1 -exp(-exp(-40.647+22.656x))
-3-2-10123
Standard Normal Deviate
-3-2-1 0 1 2 3 4
Standard Normal Deviate
Obrázek : Modely pro úmrtnost Potemníka skladištního.
Jan Koláček (PřF MU) M5VM05 Statistické modelování
21
Modely pro poissonovská data
Předpokládejme, že náhodný výběr rozsahu n je z Poissonova rozdělení, tj
Yj ... počet výskytu sledovaného jevu v určitém časovém intervalu (na pl
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).
\j > 0; y = 0,1,2,... jinak
pricemz
EYv = DYv = Ay.
velikosti t apod.).
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
Modely pro poissonovská data
Náhodnou veličinou, která má Poissonovo rozdělení, je tedy např.
• počet vadných výrobků 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ě
o 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 pol mikroskopu
• 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é modelování
Modely pro poissonovská data
Předpokládejme, že Y j závisí na k veličinách X{\,... ,Xfc, a úkolem je najít vztah mezi nimi, tj. hledáme funkci
\j = A(xz-) = \(xň,.. .,xik).
GLM modely => modelujeme pravděpodobnosti Az- pomocí linkovacích funkcí
Nejjednodušším je lineární model
A; = x?j8.
Tento model má řadu nevýhod, především je třeba zajistit, aby xjfi nabývala pouze kladných hodnot. Nejčastěji se volí tyto dvě možnosti: log-lineární model :
EYi = m = X{ = exp(x^) gl(m) = tji = gl(\i) = log(Af) = xjfi
odmocninový model :
2 /—
EYi = m = Ái = (x[j6)        =>      g2(m) = tu = g2(Ai) = V A; = x[j6
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
Modelování binomických dat pomocí poissonovského modelu
Pomocí Poissonova rozdělení Po(A) lze dobře aproximovat binomické rozdělení Bi(n,n) za podmínek
n ^ co   & 71—^0   &   nn —>> A < oo,
obvykle se doporučuje n > 30 a ti < 0,1.
Chceme-li tedy aproximovat binomické rozdělení Bi(tij, 7ij) pomocí Poissonova rozdělení Po(Az- = tijTľj) a přitom použijeme logaritmickou linkovací funkci, platí
Xi = UiTii = exp(xfjS)     log(Af) =  log(nf)  +log(7rz-) = x[j8
T
tzv. „offset"
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
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 mesic 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é modelování
26 / 33
Příklad
link function g(n)=n
	20	-1-1-1-	-.-.-.-»T-,
			•
C/)	15		•
Q			• •
<			•
			
O	10		•               • •■
			
<D		•	•
			•                                  •__■
E			• ---
	5	•	• ^^-^
		• •	
		• •	
			
		••• —	
	01		
		5           10 15	20          25          30 35
Normal Probability Plot
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
y=(0.553+0.0944xf
-3-2-101234
Standard Normal Deviate
Normal Probability Plot
<5 >
CD Q
	y
	
y	
y y	
-3-2-10123
Standard Normal Deviate Normal Probability Plot
	y y
	
y y	
y y y	
-3-2-1 0 1 2
Standard Normal Deviate
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
Úlohy k procvičení
Příklad 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é modelování
28 / 33
Úlohy k procvičení
Příklad 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é:
Infection-Severity    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.
4
Jan Koláček (PřF MU)
M5VM05 Statistické modelování
29 / 33
Úlohy k procvičení
Příklad 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ání
30 / 33
Úlohy k procvičení
Příklad 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říjmu 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é modelování
31 / 33
Úlohy k procvičení I
Příklad 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 25 m, 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ání
32 / 33
Úlohy k procvičení II
Příklad 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í p H 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é modelování
33 / 33